

###########################
### LOAD AND CLEAN DATA ###
###########################
### 1. Load packages and data ####
library(pacman)
p_load(readstata13, foreign, tidyverse, stargazer, stringr, 
       miceadds, ggplot2, reshape2, corrplot, ggpubr, dplyr, 
       coefplot, arules, sampleSelection, car,
       randomizr, estimatr, grid, gridExtra, infer, here, 
       lmtest, RItools, coin, margins, emmeans, MatchIt)


#load original data
i_am("maintext.R")
origdata<-read.dta(here("JuryDataSummer2004.dta"))

#Getting rid of non-standard jury sizes
dat<-origdata %>% filter(num_js==6) 

#Correcting jury 492
dat$ordernum<-ifelse(dat$jurynum==492,"a",dat$ordernum)
#Correcting jury 336
dat$scenario<-ifelse(dat$jurynum==336,5,dat$scenario)

#Labeling scenarios
dat$scenario<-factor(dat$scenario, labels=c("Motorcycle  1", "Circus  2", "Airbag  3", "Baldness  4", "Exercise video  5", "Childproof cap  6", "Workplace benzene  7", "Pajamas  8", "Lift-chair  9", "Computer monitor cancer  10", "Sulfur leak  11", "Train accident  12", "Escalator  13", "Parking lot abduction  14", "Trojan Yachts  15"))


#Deleting those observations w/ missing info for age, gender, race, or education 
#(losing 52 jurors):
dat<-subset.data.frame(dat,is.na(age)==F&is.na(sex)==F&is.na(white)==F&is.na(educ)==F)
#44 jurors missing income will be dealt with via dummy variable later



### 2. Remove juries with any jurors missing data ####

#Check which juries have less than 5 jurors after dropping the above jurors
toofew <- dat %>% group_by(jurynum) %>% summarize(n()) %>% filter(`n()`!=6) %>% select(jurynum)

#Drop juries with less than 6 jurors with complete info
#losing 224 jurors, 46 juries
dat <- dat %>% filter(!(jurynum %in% toofew$jurynum))
#Double check if we dropped the incomplete juries entirely.
dat %>% group_by(jurynum) %>% summarize(n()) %>% filter(`n()`!=6)

rm(toofew)

### 3. Data Cleaning, recoding, variable generation ##########

### CREATING DVs

#Scaling individual "punishment scale" decisions:
dat$scaleis<-(dat$iscale)/8
#Scaling jury-level "punishment scale" decisions:
dat$scalejr<-(dat$jryscale)/8
#Creating logged dollar amounts, adding 1 to deal with zeroes:
dat$lidoll<-log(dat$idoll + 1)
dat$ljrydoll<-log(dat$jrydoll + 1)

#create dollar amount ordinal scale 
dat <- dat %>% mutate(or_idoll = case_when(idoll==0 ~ 0,
                                           idoll%in%c(1:50000)~1,
                                           idoll%in%c(50001:100000)~2,
                                           idoll%in%c(100001:200000)~3,
                                           idoll%in%c(200001:500000-1)~4,
                                           idoll%in%c(500000:750000)~5,
                                           idoll%in%c(750001:1000000)~6,
                                           idoll%in%c(1000001:3000000)~7,
                                           idoll > 3000000 ~8,
                                           TRUE ~ NA_real_))

#rescale 0 to 1
dat$or_idoll<-dat$or_idoll/8

dat$eight.doll <- dat$or_idoll*8

#check distribution of ordinal scale dollar values
ggplot(dat, aes(x=or_idoll)) +   
  labs(title="Histogram - Ordinal Individual Dollar Punishment",x="Ordinal Scale", y = "Count")+
  geom_histogram(binwidth=0.125, color="black", fill="white") + scale_x_continuous(breaks=(0:8)/8)

## Control Variables:
#Calculating logged median pre-deliberation award, punishment, SDs:
dat <- dat %>% group_by(jurynum) %>% 
  mutate(med_or_idoll=median(or_idoll, na.rm=T),
  med_lidoll=median(lidoll, na.rm=T), 
  med_scale = median(scaleis,na.rm = T),
  sd_or_idoll = sd(or_idoll, na.rm=TRUE),
  sd_lidoll = sd(lidoll, na.rm=TRUE),
  sd_scale = sd(scaleis, na.rm=TRUE))


### Individual-level IVs used here and later:
#Female indicator:
dat$ind_female<-NA
dat$ind_female[dat$sex=="Female  2"]<-1
dat$ind_female[dat$sex=="Male  1"]<-0
#Non-white indicator:
dat$ind_nonwhite<-NA
dat$ind_nonwhite[dat$white=="Not White  0"]<-1
dat$ind_nonwhite[dat$white=="White  1"]<-0
# Hispanic indicator:
dat$ind_hisp <- ifelse(dat$ethnic=="Hispanic  4", 1, 0)
#Age indicators:
dat$ind_ageYoung<-ifelse(dat$age%in%c("18-29  1"),1,0)
dat$ind_ageMiddle<-ifelse(dat$age%in%c("30-39  2","40-49  3","50-59  4"),1,0)
dat$ind_ageOld<-ifelse(dat$age%in%c("60-69  5","70+  6"),1,0)
#Education indicators:
dat$ind_eduHS<-ifelse(dat$educ%in%c("Some high school or less  1","High school diploma  2"),1,0)
dat$ind_eduSC<-ifelse(dat$educ%in%c("Some college  3"),1,0)
dat$ind_eduCD<-ifelse(dat$educ%in%c('College diploma  4',"Some graduate school  5","Graduate degree  6"),1,0)
#Income indicators:
dat$ind_incLow<-ifelse(dat$income%in%c("$10,000 or less  1","$10,001 - $20,000  2",
                                       "$20,001 - $30,000  3"),1,0)
dat$ind_incMiddle<-ifelse(dat$income%in%c("$30,001 - $50,000  4"),1,0)
dat$ind_incHigh<-ifelse(dat$income%in%c("$50,001 - $75,000  5","over $ 75,000  6"),1,0)
dat$ind_incNA<-ifelse(is.na(dat$income)==T,1,0)


#Create jury-level medians for each gender type. 
dat <- dat %>% group_by(jurynum, ind_female) %>% mutate(gender_med_or_idoll=median(or_idoll, na.rm=T),
                                                        gender_med_scale = median(scaleis,na.rm = T))
#Create jury-level medians for each race type.
dat <- dat %>% group_by(jurynum, ind_nonwhite) %>% mutate(race_med_or_idoll =median(or_idoll, na.rm=T),
                                                          race_med_scale    =median(scaleis,na.rm = T))

#Create jury-level medians for age: young - old .
dat <- dat %>% group_by(jurynum, ind_ageOld) %>% mutate(old_med_or_idoll =median(or_idoll, na.rm=T),
                                                        old_med_scale    =median(scaleis,na.rm = T))

#Create jury-level medians for education: HS - CD .
dat <- dat %>% group_by(jurynum, ind_eduCD) %>% mutate(cd_med_or_idoll =median(or_idoll, na.rm=T),
                                                       cd_med_scale    =median(scaleis,na.rm = T))

#Create jury-level medians for income: poor - rich .
dat <- dat %>% group_by(jurynum, ind_incHigh) %>% mutate(rich_med_or_idoll =median(or_idoll, na.rm=T),
                                                         rich_med_scale    =median(scaleis,na.rm = T))

#Estimate differences between furthest-apart demographic group medians.
#missing when jury doesn't have any variation on the variable in question
dat <- dat %>% group_by(jurynum) %>% mutate(diff_race_med_or_idoll  = abs(max(race_med_or_idoll,na.rm = T)- min(race_med_or_idoll, na.rm = T)),
                                            diff_race_med_scale     = abs(max(race_med_scale ,na.rm=T)    - min(race_med_scale ,  na.rm=T)),
                                            diff_gender_med_or_idoll= abs(max(gender_med_or_idoll,na.rm=T)- min(gender_med_or_idoll ,  na.rm=T)),
                                            diff_gender_med_scale   = abs(max(gender_med_scale ,  na.rm=T)- min(gender_med_scale ,  na.rm=T)),
                                            diff_old_med_or_idoll   = abs(max(old_med_or_idoll ,  na.rm=T)- min(old_med_or_idoll ,  na.rm=T)),
                                            diff_old_med_scale      = abs(max(old_med_scale ,  na.rm=T)   - min(old_med_scale ,  na.rm=T)),
                                            diff_cd_med_or_idoll    = abs(max(cd_med_or_idoll ,  na.rm=T) - min(cd_med_or_idoll ,  na.rm=T)),
                                            diff_cd_med_scale       = abs(max(cd_med_scale ,  na.rm=T)    - min(cd_med_scale ,  na.rm=T)),
                                            diff_rich_med_or_idoll  = abs(max(rich_med_or_idoll,na.rm=T)  - min(rich_med_or_idoll ,  na.rm=T)),
                                            diff_rich_med_scale     = abs(max(rich_med_scale ,  na.rm=T)  - min(rich_med_scale ,  na.rm=T)))


###Create group-level compositional variables:
#versions including and excluding the observation's juror:

##Group-level number of women:
dat<-dat %>% group_by(jurynum) %>% mutate(totfem.inc=sum(ind_female),
                                          totfem.exc=(sum(ind_female)-ind_female))
dat$female12.inc<-ifelse(dat$totfem.inc%in%c(1,2),1,0)
dat$female3.inc<-ifelse(dat$totfem.inc%in%c(3),1,0)
dat$female4.inc<-ifelse(dat$totfem.inc%in%c(4),1,0)
dat$female56.inc<-ifelse(dat$totfem.inc%in%c(5,6),1,0)
dat$female02.exc<-ifelse(dat$totfem.exc%in%c(0,1,2),1,0)
dat$female3.exc<-ifelse(dat$totfem.exc%in%c(3),1,0)
dat$female45.exc<-ifelse(dat$totfem.exc%in%c(4,5),1,0)
prop.table(table(dat$totfem.inc))

##Group-level for number of non-white:
dat<-dat %>% group_by(jurynum) %>% mutate(totnwhite.inc=sum(ind_nonwhite),
                                          totnwhite.exc=(sum(ind_nonwhite)-ind_nonwhite))
dat$nonwhite0.inc<-ifelse(dat$totnwhite.inc%in%c(0),1,0)
dat$nonwhite1.inc<-ifelse(dat$totnwhite.inc%in%c(1),1,0)
dat$nonwhite24.inc<-ifelse(dat$totnwhite.inc%in%c(2,3,4),1,0)
dat$nonwhite0.exc<-ifelse(dat$totnwhite.exc%in%c(0),1,0)
dat$nonwhite1.exc<-ifelse(dat$totnwhite.exc%in%c(1),1,0)
dat$nonwhite24.exc<-ifelse(dat$totnwhite.exc%in%c(2,3,4),1,0)
prop.table(table(dat$totnwhite.inc))

## Group-level Hispanic:
dat <- dat %>%
  group_by(jurynum) %>%
  mutate(tothisp.inc = sum(ind_hisp),
         tothisp.exc = sum(ind_hisp) - ind_hisp)
dat$hisp0.exc <- ifelse(dat$tothisp.exc==0, 1, 0)
dat$hisp1.exc <- ifelse(dat$tothisp.exc==1, 1, 0)
dat$hisp23.exc <- ifelse(dat$tothisp.exc>1, 1, 0)
dat$hisp0.inc <- ifelse(dat$tothisp.inc==0, 1, 0)
dat$hisp1.inc <- ifelse(dat$tothisp.inc==1, 1, 0)
dat$hisp23.inc <- ifelse(dat$tothisp.inc>1, 1, 0)


##Group-level for number young:
dat<-dat %>% group_by(jurynum) %>% mutate(totyoung.inc=sum(ind_ageYoung),
                                          totyoung.exc=(sum(ind_ageYoung)-ind_ageYoung))
dat$young0.inc<-ifelse(dat$totyoung.inc%in%c(0),1,0)
dat$young1.inc<-ifelse(dat$totyoung.inc%in%c(1),1,0)
dat$young2.inc<-ifelse(dat$totyoung.inc%in%c(2),1,0)
dat$young34.inc<-ifelse(dat$totyoung.inc%in%c(3,4),1,0)
dat$young0.exc<-ifelse(dat$totyoung.exc%in%c(0),1,0)
dat$young1.exc<-ifelse(dat$totyoung.exc%in%c(1),1,0)
dat$young24.exc<-ifelse(dat$totyoung.exc%in%c(2,3,4),1,0)
prop.table(table(dat$totyoung.inc))

##Group-level for number old:
dat<-dat %>% group_by(jurynum) %>% mutate(totold.inc=sum(ind_ageOld),
                                          totold.exc=(sum(ind_ageOld)-ind_ageOld))
dat$old0.inc<-ifelse(dat$totold.inc%in%c(0),1,0)
dat$old1.inc<-ifelse(dat$totold.inc%in%c(1),1,0)
dat$old2.inc<-ifelse(dat$totold.inc%in%c(2),1,0)
dat$old35.inc<-ifelse(dat$totold.inc%in%c(3,4,5),1,0)
dat$old0.exc<-ifelse(dat$totold.exc%in%c(0),1,0)
dat$old1.exc<-ifelse(dat$totold.exc%in%c(1),1,0)
dat$old2.exc<-ifelse(dat$totold.exc%in%c(2),1,0)
dat$old35.exc<-ifelse(dat$totold.exc%in%c(3,4,5),1,0)
prop.table(table(dat$totyoung.inc))

##Group-level for number w/ hs educ or less
dat<-dat %>% group_by(jurynum) %>% mutate(toths.inc=sum(ind_eduHS),
                                          toths.exc=(sum(ind_eduHS)-(ind_eduHS)))
dat$hs0.inc<-ifelse(dat$toths.inc%in%c(0),1,0)
dat$hs1.inc<-ifelse(dat$toths.inc%in%c(1),1,0)
dat$hs2.inc<-ifelse(dat$toths.inc%in%c(2),1,0)
dat$hs35.inc<-ifelse(dat$toths.inc%in%c(3,4,5),1,0)
dat$hs0.exc<-ifelse(dat$toths.exc%in%c(0),1,0)
dat$hs1.exc<-ifelse(dat$toths.exc%in%c(1),1,0)
dat$hs2.exc<-ifelse(dat$toths.exc%in%c(2),1,0)
dat$hs35.exc<-ifelse(dat$toths.exc%in%c(3,4,5),1,0)
##Group-level for number w/ college diploma or more
dat<-dat %>% group_by(jurynum) %>% mutate(totcd.inc=sum(ind_eduCD),
                                          totcd.exc=(sum(ind_eduCD)-(ind_eduCD)))
dat$college01.inc<-ifelse(dat$totcd.inc%in%c(0,1),1,0)
dat$college2.inc<-ifelse(dat$totcd.inc%in%c(2),1,0)
dat$college3.inc<-ifelse(dat$totcd.inc%in%c(3),1,0)
dat$college45.inc<-ifelse(dat$totcd.inc%in%c(4,5),1,0)
dat$college0.exc<-ifelse(dat$totcd.exc%in%c(0),1,0)
dat$college1.exc<-ifelse(dat$totcd.exc%in%c(1),1,0)
dat$college2.exc<-ifelse(dat$totcd.exc%in%c(2),1,0)
dat$college3.exc<-ifelse(dat$totcd.exc%in%c(3),1,0)
dat$college45.exc<-ifelse(dat$totcd.exc%in%c(4,5),1,0)
##Group-level for number w/ low income
dat<-dat %>% group_by(jurynum) %>% mutate(totpoor.inc=sum(ind_incLow),
                                          totpoor.exc=(sum(ind_incLow)-(ind_incLow)))
dat$lowinc0.inc<-ifelse(dat$totpoor.inc%in%c(0),1,0)
dat$lowinc1.inc<-ifelse(dat$totpoor.inc%in%c(1),1,0)
dat$lowinc2.inc<-ifelse(dat$totpoor.inc%in%c(2),1,0)
dat$lowinc3.inc<-ifelse(dat$totpoor.inc%in%c(3),1,0)
dat$lowinc45.inc<-ifelse(dat$totpoor.inc%in%c(4,5),1,0)
dat$lowinc0.exc<-ifelse(dat$totpoor.exc%in%c(0),1,0)
dat$lowinc1.exc<-ifelse(dat$totpoor.exc%in%c(1),1,0)
dat$lowinc2.exc<-ifelse(dat$totpoor.exc%in%c(2),1,0)
dat$lowinc35.exc<-ifelse(dat$totpoor.exc%in%c(3,4,5),1,0)
##Group-level for number w/ high income
dat<-dat %>% group_by(jurynum) %>% mutate(totrich.inc=sum(ind_incHigh),
                                          totrich.exc=(sum(ind_incHigh)-(ind_incHigh)))
dat$highinc0.inc<-ifelse(dat$totrich.inc%in%c(0),1,0)
dat$highinc1.inc<-ifelse(dat$totrich.inc%in%c(1),1,0)
dat$highinc2.inc<-ifelse(dat$totrich.inc%in%c(2),1,0)
dat$highinc3.inc<-ifelse(dat$totrich.inc%in%c(3),1,0)
dat$highinc45.inc<-ifelse(dat$totrich.inc%in%c(4,5),1,0)
dat$highinc0.exc<-ifelse(dat$totrich.exc%in%c(0),1,0)
dat$highinc1.exc<-ifelse(dat$totrich.exc%in%c(1),1,0)
dat$highinc2.exc<-ifelse(dat$totrich.exc%in%c(2),1,0)
dat$highinc35.exc<-ifelse(dat$totrich.exc%in%c(3,4,5),1,0)

# DELETE IF NOT USED LATER
##Dummies for juries with roughly 50/50 opposite socio-demographic features
dat$whitenonwhite24<-ifelse(dat$totnwhite.inc%in%c(2,3,4),1,0)
dat$youngold24<-ifelse(dat$totyoung.inc%in%c(2,3,4),1,ifelse(dat$totold.inc%in%c(2,3,4),1,0))
dat$hscd24<-ifelse(dat$toths.inc%in%c(2,3,4),1,ifelse(dat$totcd.inc%in%c(2,3,4),1,0))
dat$poorrich24<-ifelse(dat$totpoor.inc%in%c(2,3,4),1,ifelse(dat$totrich.inc%in%c(2,3,4),1,0))


####


# Create jury-level ordinal dollar variable:
dat <- dat %>% mutate(jury.eight = case_when(jrydoll==0 ~ 0,
                                             jrydoll>=1&jrydoll<=50000~1,
                                             jrydoll>=50001&jrydoll<=100000~2,
                                             jrydoll>=100001&jrydoll<=200000~3,
                                             jrydoll>=200001&jrydoll<500000~4,
                                             jrydoll>=500000&jrydoll<=750000~5,
                                             jrydoll>=750001&jrydoll<=1000000~6,
                                             jrydoll>=1000001&jrydoll<=3000000~7,
                                             jrydoll > 3000000 ~8,
                                             TRUE ~ NA_real_))


##Creating new indicators of group composition
dat$white5.exc <- dat$nonwhite0.exc
dat$white4.exc <- dat$nonwhite1.exc
dat$white13.exc <- dat$nonwhite24.exc

dat$male35.exc <- dat$female02.exc
dat$male2.exc <- dat$female3.exc
dat$male01.exc <- dat$female45.exc


#Dollar deliberation first then ratings: a
#Ratings deliberation first then dollar: b

#Divide the data into pre and post deliberation rounds:
pre_del<-subset.data.frame(dat,ordernum=="b")
post_del<-subset.data.frame(dat,ordernum=="a")

##Combine ratings and dollars into one variable, by round
##Individual ratings
pre_del$scale_all_rd1<-pre_del$scaleis
pre_del$scale_all_rd2<-pre_del$or_idoll

pre_del$med_scale_rd1<-pre_del$med_scale
pre_del$med_scale_rd2<- pre_del$med_or_idoll

pre_del$sd_scale_rd1<-pre_del$sd_scale
pre_del$sd_scale_rd2<-pre_del$sd_or_idoll

post_del$scale_all_rd1<-post_del$or_idoll
post_del$scale_all_rd2<-post_del$scaleis

post_del$med_scale_rd1<-post_del$med_or_idoll
post_del$med_scale_rd2<-post_del$med_scale

post_del$sd_scale_rd1<-post_del$sd_or_idoll 
post_del$sd_scale_rd2<-post_del$sd_scale

##Jury Verdicts
pre_del$verdict_rd1 <- pre_del$scalejr
pre_del$verdict_rd2 <- pre_del$jury.eight/8

post_del$verdict_rd1 <- post_del$jury.eight/8
post_del$verdict_rd2 <- post_del$scalejr 

dat_byround<-rbind(pre_del, post_del)




#CREATE Group-level data sets:
group.data<-dplyr::select(dat_byround,c(female12.inc,female3.inc,female4.inc,female56.inc,
                                        nonwhite0.inc,nonwhite1.inc,nonwhite24.inc,
                                        hisp0.inc, hisp1.inc, hisp23.inc,
                                        young0.inc,young1.inc,young2.inc,young34.inc,old35.inc,old2.inc,old1.inc,old0.inc,
                                        hs0.inc,hs1.inc,hs2.inc,hs35.inc,college3.inc,college2.inc,college01.inc, college45.inc,
                                        lowinc0.inc,lowinc1.inc,lowinc2.inc,lowinc3.inc,lowinc45.inc,
                                        highinc45.inc,highinc3.inc,highinc2.inc,highinc1.inc,highinc0.inc,ordernum,
                                        scalejr,ljrydoll,jrydoll,med_scale,med_or_idoll,sd_scale,sd_or_idoll, med_lidoll, sd_lidoll, jurynum,scenario, 
                                        diff_race_med_or_idoll,diff_race_med_scale, diff_gender_med_or_idoll,
                                        verdict_rd1, verdict_rd2, med_scale_rd1, med_scale_rd2, sd_scale_rd1, sd_scale_rd2,
                                        diff_old_med_or_idoll, diff_old_med_scale, diff_cd_med_or_idoll, diff_cd_med_scale, diff_rich_med_or_idoll, diff_rich_med_scale,
                                        diff_gender_med_scale, totfem.inc, totnwhite.inc, totold.inc, totcd.inc, totrich.inc,totpoor.inc, toths.inc, totyoung.inc, jury.eight))
group.data<-unique(group.data)
group.data$or_jrydoll<- group.data$jury.eight/8

#Separate dollar-first jury data from punishment scale-first jury data.
group.pre_del<-subset.data.frame(group.data,ordernum=="b")
group.post_del<-subset.data.frame(group.data,ordernum=="a")

#check scenario mean punishment scales
post_del<-post_del %>% group_by(scenario) %>% mutate(meanscaleis=mean(scaleis, na.rm=TRUE))
scenariomeans<- unique(post_del[,c("meanscaleis", "scenario")]) 
as.data.frame(scenariomeans[order(scenariomeans$meanscaleis),])

#Make Workplace benzene  7 the base category in factor variable "scenario".
dat_byround <- dat_byround %>%
  mutate(scenario = relevel(scenario, ref = "Workplace benzene  7"))

post_del <- post_del %>%
  mutate(scenario = relevel(scenario, ref = "Workplace benzene  7"))
pre_del <- pre_del %>%
  mutate(scenario = relevel(scenario, ref = "Workplace benzene  7"))


group.post_del <- group.post_del %>%
  mutate(scenario = relevel(scenario, ref = "Workplace benzene  7"))
group.pre_del <- group.pre_del %>%
  mutate(scenario = relevel(scenario, ref = "Workplace benzene  7"))

group.data <- group.data %>%
  mutate(scenario = relevel(scenario, ref = "Workplace benzene  7"))



#Create hung jury dependent variables: nothung=1, hung=0
group.post_del$or_jrydoll_nothung<-1
group.post_del$or_jrydoll_nothung[is.na(group.post_del$jury.eight)==TRUE]<-0

group.post_del$scalejr_nothung<-1
group.post_del$scalejr_nothung[is.na(group.post_del$scalejr)==TRUE]<-0

group.pre_del$or_jrydoll_nothung<-1
group.pre_del$or_jrydoll_nothung[is.na(group.pre_del$jury.eight)==TRUE]<-0

group.pre_del$scalejr_nothung<-1
group.pre_del$scalejr_nothung[is.na(group.pre_del$scalejr)==TRUE]<-0

group.pre_del$round1_nothung<-1
group.pre_del$round1_nothung[is.na(group.pre_del$scalejr)==TRUE]<-0

group.pre_del$round2_nothung<-1
group.pre_del$round2_nothung[is.na(group.pre_del$jury.eight)==TRUE]<-0

group.post_del$round1_nothung<-1
group.post_del$round1_nothung[is.na(group.post_del$jury.eight)==TRUE]<-0

group.post_del$round2_nothung<-1
group.post_del$round2_nothung[is.na(group.post_del$scalejr)==TRUE]<-0

# race difference medians

group.pre_del$diff_race_med_round1 <- group.pre_del$diff_race_med_scale
group.pre_del$diff_race_med_round2 <- group.pre_del$diff_race_med_or_idoll

group.post_del$diff_race_med_round1 <- group.post_del$diff_race_med_or_idoll
group.post_del$diff_race_med_round2 <- group.post_del$diff_race_med_scale



group.data$or_jrydoll_nothung<-1
group.data$or_jrydoll_nothung[is.na(group.data$jury.eight)==TRUE]<-0

group.data$scalejr_nothung<-1
group.data$scalejr_nothung[is.na(group.data$scalejr)==TRUE]<-0



percentiles<- c(1, 50000, 100000, 200000, 500000-1, 750000, 1000000, 3000000)
group.percentiles<- c(1, 50000, 100000, 200000, 500000-1, 750000, 1000000, 3000000)

group.dat_byround<-rbind(group.pre_del, group.post_del)

toofew_byround <- dat_byround %>% group_by(jurynum) %>% summarize(n()) %>% filter(`n()`!=6) %>% select(jurynum)
rm(toofew_byround)
rm(percentiles, group.percentiles, scenariomeans)


###########################
### ANALYSIS: MAIN TEXT ###
###########################



## Table 1. Predictors of Pre-Deliberation Individual Punitiveness #####

# Column 1: Individual Predictors Only
predelib_prefs_ind<-lm(scale_all_rd1~
                         ind_nonwhite+
                         ind_female+
                         ind_ageYoung+ind_ageOld+
                         ind_eduHS+ind_eduCD+
                         ind_incLow+ind_incHigh+
                         ind_incNA+
                         scenario,data=dat_byround)
# Column 2: add composition variables
predelib_prefs_grp<-lm(scale_all_rd1~
                         ind_nonwhite+nonwhite1.exc+nonwhite0.exc+
                         ind_female+female3.exc+female02.exc+
                         ind_ageYoung+ind_ageOld+old1.exc+old2.exc+old35.exc+
                         ind_eduHS+ind_eduCD+college1.exc+college2.exc+college3.exc+college45.exc+
                         ind_incLow+ind_incHigh+highinc1.exc+highinc2.exc+highinc35.exc+
                         ind_incNA+
                         scenario,data=dat_byround)

var_order <- c("ind_nonwhite", "ind_female", "ind_ageYoung", "ind_ageOld", "ind_eduHS", "ind_eduCD", "ind_incLow", "ind_incHigh", "ind_incNA", "nonwhite1.exc", "nonwhite0.exc", "female3.exc", "female02.exc", "old1.exc", "old2.exc", "old35.exc", "college1.exc", "college2.exc", "college3.exc", "college45.exc", "highinc1.exc", "highinc2.exc", "highinc35.exc")

stargazer( predelib_prefs_ind,  predelib_prefs_grp, 
           se=list(coef(summary(predelib_prefs_ind,cluster = c("jurynum")))[, 2],
                   coef(summary(predelib_prefs_grp,cluster = c("jurynum")))[, 2]),
           type="latex", out="Final Models/predelib_prefs_all.tex", star.char = c("*", "**", "***"), star.cutoffs = c(.05, .01, .001)
           , omit= "scenario"
           , order=(var_order)
           ,covariate.labels = c("POC" ,"Female" ,  "Young" ,  "Older" , 
                                 "High School Grad" ,  "College Grad" , 
                                 "Low Income" ,  "High Income" , "Missing Income" ,
                                 "4 Whites", "5 Whites", 
                                 "2 Men", "3+ Men", 
                                 "1 Older", "2 Older", "3+ Older",  
                                 "1 College Grad", "2 College Grad", "3 College Grad", "4+ College Grad",   
                                 "1 High Income", "2 High Income", "3+ High Income")
           , omit.stat=c("f", "ser"), omit.table.layout = "n", dep.var.labels=c("Pre-Deliberation Preferences (t_0)"), no.space = T, single.row = T,
           add.lines=list(c("Legal Case Fixed Effects?", "Yes", "Yes")))


##Footnote 18: F-test on models with only group variables
predelib_prefs_grp_only<-lm(scale_all_rd1~
                              nonwhite1.exc+nonwhite0.exc+
                              female3.exc+female02.exc+
                              old1.exc+old2.exc+old35.exc+
                              college1.exc+college2.exc+college3.exc+college45.exc+
                              highinc1.exc+highinc2.exc+highinc35.exc,
                            data=dat_byround)
summary(predelib_prefs_grp_only)

rm(predelib_prefs_grp, predelib_prefs_grp_only, predelib_prefs_ind)

## Table 2. Predictors of Post-Deliberation Individual Punitiveness #####

## Column 1: Race variables only
postdelib_prefs_simp <- lm(scale_all_rd2~
                             ind_nonwhite+nonwhite0.exc+nonwhite1.exc+
                             scale_all_rd1+
                             scenario,data=dat_byround)

## Column 2: All group-level variables + pre-deliberation control
postdelib_prefs_all<-lm(scale_all_rd2~
                          ind_nonwhite+nonwhite0.exc+nonwhite1.exc+
                          ind_female+female02.exc+female3.exc+
                          ind_ageYoung+ind_ageOld+old35.exc+old2.exc+old1.exc+
                          ind_eduHS+ind_eduCD+college45.exc+college3.exc+college2.exc+college1.exc+
                          ind_incLow+ind_incHigh+highinc35.exc+highinc2.exc+highinc1.exc+
                          ind_incNA+scale_all_rd1+
                          scenario,data=dat_byround)

## Column 3: all group-level variables + pre-deliberation control + pre-deliberation group median
postdelib_prefs_all_2<-lm(scale_all_rd2~
                            ind_nonwhite+nonwhite0.exc+nonwhite1.exc+
                            ind_female+female02.exc+female3.exc+
                            ind_ageYoung+ind_ageOld+old35.exc+old2.exc+old1.exc+
                            ind_eduHS+ind_eduCD+college45.exc+college3.exc+college2.exc+college1.exc+
                            ind_incLow+ind_incHigh+highinc35.exc+highinc2.exc+highinc1.exc+
                            ind_incNA+scale_all_rd1+med_scale_rd1+
                            scenario,data=dat_byround)

##Create table
stargazer(postdelib_prefs_simp, postdelib_prefs_all, postdelib_prefs_all_2, 
          se=list(coef(summary(postdelib_prefs_simp,cluster = c("jurynum")))[, 2],
                  coef(summary(postdelib_prefs_all,cluster = c("jurynum")))[, 2],
                  coef(summary(postdelib_prefs_all_2,cluster = c("jurynum")))[, 2]),
          type="latex", out="Final Models/postdelib_main_all.tex", star.char = c("*", "**", "***"), star.cutoffs = c(.05, .01, .001)
          , omit= c("scenario", "ind_nonwhite", "ind_female", 
                    "ind_ageYoung", "ind_ageOld", "ind_eduHS", "ind_eduCD",
                    "ind_incLow", "ind_incHigh", "ind_incNA")
          , order=(var_order)
          ,covariate.labels = c("4 Whites", "5 Whites", 
                                "2 Men", "3+ Men", 
                                "1 Older", "2 Older", "3+ Older",  
                                "1 College Grad", "2 College Grad", "3 College Grad", "4+ College Grad",   
                                "1 High Income", "2 High Income", "3+ High Income",
                                "Individual Pre-delib. Preference" ,  "Jury Median Pre-delib. Preference")
          , omit.stat=c("f", "ser"), dep.var.labels=c("Post-Deliberation Preference (t2)"),
          add.lines = list(c("Legal Case Fixed Effects?", "Yes", "Yes", "Yes"),
                           c("Indiv. Demographic Controls?", "Yes", "Yes", "Yes")),
          column.sep.width = "2pt", omit.table.layout = "n", no.space = T)


rm(postdelib_prefs_all_2, postdelib_prefs_simp)

##Footnote 19: Correlation between jury's median predeliberation preferences and the number of nonwhites in the group 
cor(group.dat_byround$med_scale_rd1, group.dat_byround$totnwhite.inc)

####


## Figure 2. Predicted Individual Punitiveness by Individual and Group Race #####

#Interact composition with race for figure
postdelib_interact_race_all<-lm(scale_all_rd2~
                                  ind_nonwhite*nonwhite1.exc+ind_nonwhite*nonwhite0.exc+
                                  ind_female+female3.exc+female02.exc+
                                  ind_ageYoung+ind_ageOld+old1.exc+old2.exc+old35.exc+
                                  ind_eduHS+ind_eduCD+college1.exc+college2.exc+college3.exc+college45.exc+
                                  ind_incLow+ind_incHigh+highinc1.exc+highinc2.exc+highinc35.exc+
                                  ind_incNA+scale_all_rd1+scenario,data=dat_byround)

##### Predicted scores based on the base models: ratings + dollars

# create toy data frame for predictions
pp_all.data <- data.frame(
  "ind_nonwhite"=c(1,0,1,0,1,0),
  "nonwhite0.exc"=c(1,1,0,0,0,0),
  "nonwhite1.exc"=c(0,0,1,1,0,0),
  "ind_female"=rep(median(dat_byround$ind_female),6),
  "female02.exc"=rep(0,6),
  "female3.exc"=rep(1,6),
  "ind_ageOld"= rep(median(dat_byround$ind_ageOld),6),
  "ind_ageYoung"= rep(median(dat_byround$ind_ageYoung),6),
  "old35.exc"= rep(0,6),
  "old1.exc"= rep(1,6), 
  "old2.exc"= rep(0,6),
  "ind_eduCD"= rep(median(dat_byround$ind_eduCD),6),
  "ind_eduHS"= rep(median(dat_byround$ind_eduHS),6),
  "college45.exc"= rep(0,6),
  "college1.exc"= rep(0,6),
  "college2.exc"= rep(1,6),
  "college3.exc"= rep(0,6),
  "ind_incHigh"= rep(median(dat_byround$ind_incHigh),6),
  "ind_incLow"= rep(median(dat_byround$ind_incLow),6),
  "highinc35.exc" =rep(0,6),
  "highinc1.exc" = rep(0,6),
  "highinc2.exc" = rep(1,6),
  "ind_incNA"= rep(0,6),
  "scale_all_rd1"=rep(median(dat_byround$scale_all_rd1[is.na(dat_byround$scale_all_rd1)==F]),6),
  "scenario"= rep("Workplace benzene  7",6)  
)

# calculate predicted values
pp_all.values<-as.data.frame(predict(postdelib_interact_race_all,pp_all.data,se.fit = T))

# add confidence intervals 
pp_all.values <- mutate(pp_all.values,
                        lwr=fit-1.39*se.fit,
                        upr=fit+1.39*se.fit)

# combine data for graph
pp_all.graph.data<-data.frame("values"=pp_all.values$fit,
                              "ind_nonwhite"=c("POC","White","POC","White","POC","White"),
                              "group_nonwhite.exc"=c("0","0","1","1","2-4","2-4"),
                              "lwr"=pp_all.values$lwr,
                              "upr"=pp_all.values$upr)
pp_all.graph.data$race<-factor(pp_all.graph.data$ind_nonwhite,levels=c("POC","White"))

#recode for graph labeling
pp_all.graph.data$group_nonwhite.exc<-as.factor(pp_all.graph.data$group_nonwhite.exc)
pp_all.graph.data$group_nonwhite.exc<-factor(pp_all.graph.data$group_nonwhite.exc, levels=c("2-4", "1", "0"))
pp_all.graph.data$ind_nonwhite <- factor(pp_all.graph.data$ind_nonwhite, labels=c("POC Jurors", "White Jurors"))

#produce plot:
nonwhiteplot_all_facet<-ggplot(pp_all.graph.data, aes(x = group_nonwhite.exc, y = values, fill = ind_nonwhite))+ 
  geom_bar(stat="identity", width=.5, position = "dodge")+
  geom_pointrange(aes(x=group_nonwhite.exc,ymin=lwr,ymax=upr,group=ind_nonwhite),colour="black",
                  position=position_dodge(width = .5))+
  labs(x="Number of White Jurors in Group, Excluding Respondent",y="Predicted Post-Delib. Preference")+
  theme_bw()+
  scale_fill_grey(start = .4, end = .8,name="Juror Race")+
  scale_x_discrete(labels=c("3 or Fewer","4","5"))+
  facet_wrap(~ind_nonwhite) +
  theme(
    legend.position="none",
    rect = element_rect(fill = "transparent"),
    axis.text.x = element_text(colour = "black",face="bold"),
    axis.text.y = element_text(colour = "black",face="bold")
  )+
  ylim(0, 0.9)
nonwhiteplot_all_facet
ggsave(here("Figures/predicted_punishment_all_facet.pdf"), width = 10, height = 8)

rm(nonwhiteplot_all_facet, pp_all.data, pp_all.graph.data, pp_all.values)

## Table 3: Predicting hung juries and verdicts #####

## Run selection models:
  #round 1
heckit_withPM_controls_rd1_all<-heckit( round1_nothung~
                                          nonwhite0.inc+nonwhite1.inc+
                                          female12.inc+female3.inc+female4.inc+
                                          old35.inc+old2.inc+old1.inc+
                                          college45.inc+college3.inc+college2.inc+
                                          highinc45.inc+highinc3.inc+highinc2.inc+highinc1.inc+
                                          scenario+sd_scale_rd1,
                                        
                                        verdict_rd1~
                                          nonwhite0.inc+nonwhite1.inc+
                                          female12.inc+female3.inc+female4.inc+
                                          old35.inc+old2.inc+old1.inc+
                                          college45.inc+college3.inc+college2.inc+
                                          highinc45.inc+highinc3.inc+highinc2.inc+highinc1.inc+
                                          scenario+med_scale_rd1,data=group.dat_byround)
  #round 2
heckit_withPM_controls_rd2_all<-heckit( round2_nothung~
                                          nonwhite0.inc+nonwhite1.inc+
                                          female12.inc+female3.inc+female4.inc+
                                          old35.inc+old2.inc+old1.inc+
                                          college45.inc+college3.inc+college2.inc+
                                          highinc45.inc+highinc3.inc+highinc2.inc+highinc1.inc+
                                          scenario+sd_scale_rd1,
                                        
                                        verdict_rd2~
                                          nonwhite0.inc+nonwhite1.inc+
                                          female12.inc+female3.inc+female4.inc+
                                          old35.inc+old2.inc+old1.inc+
                                          college45.inc+college3.inc+college2.inc+
                                          highinc45.inc+highinc3.inc+highinc2.inc+highinc1.inc+
                                          scenario+med_scale_rd1,data=group.dat_byround)


var_order_grp <- c("nonwhite1.inc", "nonwhite0.inc", "female4.inc", "female3.inc", "female12.inc", "old1.inc", "old2.inc", "old35.inc", "college45.inc", "college3.inc", "college2.inc", "college45.inc", "highinc1.inc", "highinc2.incc", "highinc3.inc", "highinc45.inc")

# export first stage models to go in appendix
stargazer(heckit_withPM_controls_rd1_all,heckit_withPM_controls_rd2_all,  
type="latex", selection.equation=TRUE, out="Appendix Figures/heckit_st1_withPM_all.tex", star.char = c("*", "**", "***"), star.cutoffs = c(.05, .01, .001)
, omit= "scenario"
, order=(var_order_grp)
,covariate.labels = c("5 White", "6 White", 
                      "2 Men", "3 Men", "4+ Men",
                      "1 Older", "2 Older", "3+ Older",
                      "2 College Grad", "3 College Grad", "4+ College Grad", 
                      "1 High Income", "2 High Income", "3 High Income", "4+ High Income",
                      "SD Pre-Delib. Preference (t0)")
, omit.stat=c("f", "ser"), dep.var.labels=c("Reached Verdict (Round 1)", "Reached Verdict (Round 2)"),
single.row=F, column.sep.width = "2pt", omit.table.layout = "n", no.space = T
) 

#second stage models to go in main text 
stargazer(heckit_withPM_controls_rd1_all,heckit_withPM_controls_rd2_all, 
          type="latex", selection.equation=FALSE, out="Final Models/heckit_st2_withPM_all.tex", star.char = c("*", "**", "***"), star.cutoffs = c(.05, .01, .001)
          , omit= "scenario"
          , order=(var_order_grp)
          ,covariate.labels = c("5 White", "6 White", 
                                "2 Men", "3 Men", "4+ Men",
                                "1 Older", "2 Older", "3+ Older",
                                "2 College Grad", "3 College Grad", "4+ College Grad", 
                                "1 High Income", "2 High Income", "3 High Income", "4+ High Income", 
                                "Median Pre-Delib. Preference (t0)")
          , omit.stat=c("f", "ser"), dep.var.labels=c("Verdict (Round 1)", "Verdict (Round 2)")  ,
          single.row = F, column.sep.width = "2pt", omit.table.layout = "n", no.space = T
          
)

rm(heckit_withPM_controls_rd1_all, heckit_withPM_controls_rd2_all)

## Setup for figures 3 and 4: one dissenter analyses #####
  ## Create datasets for each group: actual and placebo race ####

# Subset to 1-nonwhite juries:
pre_del.hetr <- pre_del %>% filter(totnwhite.inc==1)
post_del.hetr <- post_del %>% filter(totnwhite.inc==1)

#grab all-white juries
pre_del.homog <- pre_del %>% filter(totnwhite.inc==0)
post_del.homog <- post_del %>% filter(totnwhite.inc==0)

#create df of between-race differences for each order/scale:

#dollar deliberations, dollar first juries:
# calculate differences between race groups
dr1 <- post_del.hetr %>% 
  group_by(jurynum,ind_nonwhite,jury.eight,scenario) %>% #group by jury, race; keep jury dollar decision, scenario
  summarize(mean = mean(eight.doll,na.rm=T)) #find mean of dollar prefs
#grab only whites
dr1.white <- dr1 %>% filter(ind_nonwhite==0)
#grab only nonwhites
dr1.nonwhite <- dr1 %>% filter(ind_nonwhite==1)
#merge whites and nonwhites by jury
#create variable for differences between white average and nonwhite
dr1 <- left_join(dr1.white,dr1.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
rm(dr1.nonwhite, dr1.white)
#find absolute difference
dr1 <- dr1 %>% mutate(diff=abs(diff))
#recode NaNs to NAs
dr1$diff[is.nan(dr1$diff)]<-NA
dr1$mean.x[is.nan(dr1$mean.x)]<-NA
dr1$mean.y[is.nan(dr1$mean.y)]<-NA

# repeat for dollar delib, rating first groups:
dr2 <- pre_del.hetr %>% group_by(jurynum,ind_nonwhite,jury.eight,scenario) %>% 
  summarize(mean = mean(eight.doll,na.rm=T))
dr2.white <- dr2 %>% filter(ind_nonwhite==0)
dr2.nonwhite <- dr2 %>% filter(ind_nonwhite==1)
dr2 <- left_join(dr2.white,dr2.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
rm(dr2.nonwhite, dr2.white)

dr2 <- dr2 %>% mutate(diff=abs(diff))
dr2$diff[is.nan(dr2$diff)]<-NA
dr2$mean.x[is.nan(dr2$mean.x)]<-NA
dr2$mean.y[is.nan(dr2$mean.y)]<-NA

# repeat for ratings, ratings first
pr1 <- pre_del.hetr %>% group_by(jurynum,ind_nonwhite,jryscale,scenario) %>% 
  summarize(mean = mean(iscale,na.rm=T))
pr1.white <- pr1 %>% filter(ind_nonwhite==0)
pr1.nonwhite <- pr1 %>% filter(ind_nonwhite==1)
pr1 <- left_join(pr1.white,pr1.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
rm(pr1.white, pr1.nonwhite)

pr1 <- pr1 %>% mutate(diff=abs(diff))
pr1$diff[is.nan(pr1$diff)]<-NA
pr1$mean.x[is.nan(pr1$mean.x)]<-NA
pr1$mean.y[is.nan(pr1$mean.y)]<-NA

#repeat for ratings, dollars first
pr2 <- post_del.hetr %>% group_by(jurynum,ind_nonwhite,jryscale,scenario) %>% 
  summarize(mean = mean(iscale,na.rm=T))
pr2.white <- pr2 %>% filter(ind_nonwhite==0)
pr2.nonwhite <- pr2 %>% filter(ind_nonwhite==1)
pr2 <- left_join(pr2.white,pr2.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
rm(pr2.nonwhite, pr2.white)

pr2 <- pr2 %>% mutate(diff=abs(diff))
pr2$diff[is.nan(pr2$diff)]<-NA
pr2$mean.x[is.nan(pr2$mean.x)]<-NA
pr2$mean.y[is.nan(pr2$mean.y)]<-NA

##Mean differences for each group
summary(dr1$diff)
summary(dr2$diff)
summary(pr1$diff)
summary(pr2$diff)

# create dataframes for each group

dr1.a <- dr1 %>%
  #restrict to juries where white-nonwhite diff at least 1.5
  filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
  #calculate diff from jury mean for whites and nonwhites
  mutate(diff.white = abs(mean.x-jury.eight.x),
         diff.nonwhite = abs(mean.y-jury.eight.x),
         #calculate diff between whites and nonwhites' differences from mean
         diff.diff = diff.nonwhite-diff.white)
dr2.a <- dr2 %>% 
  filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jury.eight.x),
         diff.nonwhite = abs(mean.y-jury.eight.x),
         diff.diff = diff.nonwhite-diff.white)

pr1.a <- pr1 %>% 
  filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jryscale.x),
         diff.nonwhite = abs(mean.y-jryscale.x),
         diff.diff = diff.nonwhite-diff.white)

pr2.a <- pr2 %>% 
  filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jryscale.x),
         diff.nonwhite = abs(mean.y-jryscale.x),
         diff.diff = diff.nonwhite-diff.white)

#combine and save out 
dr1.a$round <- 1
dr1.a$type <- "dollar"
dr2.a$round <- 2
dr2.a$type <- "dollar"
pr1.a$round <- 1
pr1.a$type <- "rating"
pr2.a$round <- 2
pr2.a$type <- "rating"

nonwhite_rep_actual <- rbind(dr1.a, dr2.a, pr1.a, pr2.a)
nonwhite_rep_actual$nonwhite <- 1

set.seed(9904)

#within juries, randomly assign one nonwhite juror
post_del.homog <- within(post_del.homog, {
  pseudo_id <- block_ra(blocks = post_del.homog$jurynum, prob_each = c(.83, .17))
})
pre_del.homog <- within(pre_del.homog, {
  pseudo_id <- block_ra(blocks = pre_del.homog$jurynum, prob_each = c(.83, .17))
})


##Calculate raw differences, then abs. val for each group

#Dollar Scale; Dollar-first (rd 1):
dr1_placebo <- post_del.homog %>% group_by(jurynum,pseudo_id,jury.eight,scenario) %>% 
  summarize(mean = mean(eight.doll,na.rm=T))
dr1_placebo.white <- dr1_placebo %>% filter(pseudo_id==0)
dr1_placebo.nonwhite <- dr1_placebo %>% filter(pseudo_id==1)
dr1_placebo <- left_join(dr1_placebo.white,dr1_placebo.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
dr1_placebo <- dr1_placebo %>% mutate(diff=abs(diff))

dr1_placebo.a <- dr1_placebo %>%
  filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jury.eight.x),
         diff.nonwhite = abs(mean.y-jury.eight.x),
         diff.diff = diff.nonwhite-diff.white)

#Dollar Scale; Rating-first (rd 2):
dr2_placebo <- pre_del.homog %>% group_by(jurynum,pseudo_id,jury.eight,scenario) %>% 
  summarize(mean = mean(eight.doll,na.rm=T))
dr2_placebo.white <- dr2_placebo %>% filter(pseudo_id==0)
dr2_placebo.nonwhite <- dr2_placebo %>% filter(pseudo_id==1)
dr2_placebo <- left_join(dr2_placebo.white,dr2_placebo.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
dr2_placebo <- dr2_placebo %>% mutate(diff=abs(diff))

dr2_placebo.a <- dr2_placebo %>% 
  filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jury.eight.x),
         diff.nonwhite = abs(mean.y-jury.eight.x),
         diff.diff = diff.nonwhite-diff.white)

#Rating Scale; Rating-first (rd 1):
pr1_placebo <- pre_del.homog %>% group_by(jurynum,pseudo_id,jryscale,scenario) %>% 
  summarize(mean = mean(iscale,na.rm=T))
pr1_placebo.white <- pr1_placebo %>% filter(pseudo_id==0)
pr1_placebo.nonwhite <- pr1_placebo %>% filter(pseudo_id==1)
pr1_placebo <- left_join(pr1_placebo.white,pr1_placebo.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
pr1_placebo <- pr1_placebo %>% mutate(diff=abs(diff))

pr1_placebo.a <- pr1_placebo %>% 
  filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jryscale.x),
         diff.nonwhite = abs(mean.y-jryscale.x),
         diff.diff = diff.nonwhite-diff.white)

#Rating Scale; Dollar-first (rd 2):
pr2_placebo <- post_del.homog %>% group_by(jurynum,pseudo_id,jryscale,scenario) %>% 
  summarize(mean = mean(iscale,na.rm=T))
pr2_placebo.white <- pr2_placebo %>% filter(pseudo_id==0)
pr2_placebo.nonwhite <- pr2_placebo %>% filter(pseudo_id==1)
pr2_placebo <- left_join(pr2_placebo.white,pr2_placebo.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
pr2_placebo <- pr2_placebo %>% mutate(diff=abs(diff))

pr2_placebo.a <- pr2_placebo %>% 
  filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jryscale.x),
         diff.nonwhite = abs(mean.y-jryscale.x),
         diff.diff = diff.nonwhite-diff.white)


placebo_rep <- rbind(dr1_placebo.a, dr2_placebo.a, pr1_placebo.a, pr2_placebo.a)
rm(dr1_placebo.white, dr2_placebo.white, pr1_placebo.white, pr2_placebo.white,
   dr1_placebo.nonwhite, dr2_placebo.nonwhite, pr1_placebo.nonwhite, pr2_placebo.nonwhite)


## Figure 4: Estimated Influence Gaps for Lone POC and Placebo White Dissenters ####
##Identify all dissenters in all-white groups
placebo_full <- data.frame(matrix(ncol = 16, nrow =0))

#randomly label each juror in a jury 1-6
placebo_full_rate1<-pre_del.homog %>% group_by(jurynum) %>% 
  mutate(juror_num=sample(1:6, 6, replace=F))
placebo_full_doll1<-post_del.homog %>% group_by(jurynum) %>% 
  mutate(juror_num=sample(1:6, 6, replace=F))

# treat each juror in turn as the randomly-selected placebo nonwhite juror
j <- 1
for (j in 1:6) {
  
  #Dollar Scale; Dollar-first (rd 1):
    #label the selected juror
  placebo_full_doll1$temp_rand_juror <- 0
  placebo_full_doll1$temp_rand_juror[placebo_full_doll1$juror_num==j] <- 1
    #calculate diff between placebo juror and others
  dr1_placebo_full <- placebo_full_doll1 %>% group_by(jurynum,temp_rand_juror,jury.eight,scenario) %>% 
    summarize(mean = mean(eight.doll,na.rm=T))
    #separate sets for placebo and other jurors
  dr1_placebo.white_full <- dr1_placebo_full %>% filter(temp_rand_juror==0)
  dr1_placebo.nonwhite_full <- dr1_placebo_full %>% filter(temp_rand_juror==1)
    #join together and calculate difference
  dr1_placebo_full <- left_join(dr1_placebo.white_full,dr1_placebo.nonwhite_full,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  dr1_placebo_full <- dr1_placebo_full %>% mutate(diff=abs(diff))
    #replace NaNs with NAs
  dr1_placebo_full$diff[is.nan(dr1_placebo_full$diff)]<-NA
  dr1_placebo_full$mean.x[is.nan(dr1_placebo_full$mean.x)]<-NA
  dr1_placebo_full$mean.y[is.nan(dr1_placebo_full$mean.y)]<-NA
  
  # add contextual info 
  dr1_placebo_full$loop <- j
  dr1_placebo_full$round <-1
  dr1_placebo_full$type <- "dollar"
  dr1_placebo_full$num_disagree <- 0
  dr1_placebo_full$num_disagree[dr1_placebo_full$diff>=1.5] <- 1
  
  #repeat for dollar scale; rating-first (rd 2):
  placebo_full_rate1$temp_rand_juror <- 0
  placebo_full_rate1$temp_rand_juror[placebo_full_doll1$juror_num==j] <- 1
  dr2_placebo_full <- placebo_full_rate1 %>% group_by(jurynum,temp_rand_juror,jury.eight,scenario) %>% 
    summarize(mean = mean(eight.doll,na.rm=T))
  dr2_placebo.white_full <- dr2_placebo_full %>% filter(temp_rand_juror==0)
  dr2_placebo.nonwhite_full <- dr2_placebo_full %>% filter(temp_rand_juror==1)
  dr2_placebo_full <- left_join(dr2_placebo.white_full,dr2_placebo.nonwhite_full,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  dr2_placebo_full <- dr2_placebo_full %>% mutate(diff=abs(diff))
  
  dr2_placebo_full$loop <- j
  dr2_placebo_full$round <- 2
  dr2_placebo_full$type <- "dollar"
  dr2_placebo_full$num_disagree <- 0
  dr2_placebo_full$num_disagree[dr2_placebo_full$diff>=1.5] <- 1
  table(dr2_placebo_full$num_disagree)  
  
  #repeat for rating Scale; rating-first (rd 1):
  pr1_placebo_full <- placebo_full_rate1 %>% group_by(jurynum,temp_rand_juror,jryscale,scenario) %>% 
    summarize(mean = mean(iscale,na.rm=T))
  pr1_placebo.white_full <- pr1_placebo_full %>% filter(temp_rand_juror==0)
  pr1_placebo.nonwhite_full <- pr1_placebo_full %>% filter(temp_rand_juror==1)
  pr1_placebo_full <- left_join(pr1_placebo.white_full,pr1_placebo.nonwhite_full,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  pr1_placebo_full <- pr1_placebo_full %>% mutate(diff=abs(diff))
  
  pr1_placebo_full$loop <- j
  pr1_placebo_full$round <- 1
  pr1_placebo_full$type <- "rating"
  pr1_placebo_full$num_disagree <- 0
  pr1_placebo_full$num_disagree[pr1_placebo_full$diff>=1.5] <- 1
  table(pr1_placebo_full$num_disagree)  
  
  
  #Rating Scale; Dollar-first (rd 2):
  pr2_placebo_full <- placebo_full_doll1 %>% group_by(jurynum,temp_rand_juror,jryscale,scenario) %>% 
    summarize(mean = mean(iscale,na.rm=T))
  pr2_placebo.white_full <- pr2_placebo_full %>% filter(temp_rand_juror==0)
  pr2_placebo.nonwhite_full <- pr2_placebo_full %>% filter(temp_rand_juror==1)
  pr2_placebo_full <- left_join(pr2_placebo.white_full,pr2_placebo.nonwhite_full,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  pr2_placebo_full <- pr2_placebo_full %>% mutate(diff=abs(diff))
  
  pr2_placebo_full$loop <- j
  pr2_placebo_full$round <- 2
  pr2_placebo_full$type <- "rating"
  pr2_placebo_full$num_disagree <- 0
  pr2_placebo_full$num_disagree[pr2_placebo_full$diff>=1.5] <- 1
  table(pr2_placebo_full$num_disagree)
  
  ##Create dissenters-only dataset
  placebo_full_dr1_temp <- dr1_placebo_full %>% 
    filter(num_disagree==1 & is.na(jury.eight.x)==F) %>%  # keep only dissenters 
    mutate(diff.white = abs(mean.x-jury.eight.x), #calculate differences
           diff.nonwhite = abs(mean.y-jury.eight.x),
           diff.diff = diff.nonwhite-diff.white)
  
  placebo_full_dr2_temp <- dr2_placebo_full %>% 
    filter(num_disagree==1 & is.na(jury.eight.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jury.eight.x),
           diff.nonwhite = abs(mean.y-jury.eight.x),
           diff.diff = diff.nonwhite-diff.white)
  
  placebo_full_pr1_temp <- pr1_placebo_full %>% 
    filter(num_disagree==1 & is.na(jryscale.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jryscale.x),
           diff.nonwhite = abs(mean.y-jryscale.x),
           diff.diff = diff.nonwhite-diff.white)
  
  placebo_full_pr2_temp <- pr2_placebo_full %>% 
    filter(num_disagree==1 & is.na(jryscale.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jryscale.x),
           diff.nonwhite = abs(mean.y-jryscale.x),
           diff.diff = diff.nonwhite-diff.white)
  
  
  if(j==1){
    placebo_full_1 <- rbind(placebo_full_dr1_temp, placebo_full_dr2_temp, placebo_full_pr1_temp, placebo_full_pr2_temp)
  }
  
  if (j==2){
    placebo_full_2 <- rbind(placebo_full_dr1_temp, placebo_full_dr2_temp, placebo_full_pr1_temp, placebo_full_pr2_temp)
  }
  
  if (j==3){
    placebo_full_3 <- rbind(placebo_full_dr1_temp, placebo_full_dr2_temp, placebo_full_pr1_temp, placebo_full_pr2_temp)
  }
  
  if (j==4){
    placebo_full_4 <- rbind(placebo_full_dr1_temp, placebo_full_dr2_temp, placebo_full_pr1_temp, placebo_full_pr2_temp)
  }
  
  if (j==5){
    placebo_full_5 <- rbind(placebo_full_dr1_temp, placebo_full_dr2_temp, placebo_full_pr1_temp, placebo_full_pr2_temp)
  }
  
  if (j==6){
    placebo_full_6 <- rbind(placebo_full_dr1_temp, placebo_full_dr2_temp, placebo_full_pr1_temp, placebo_full_pr2_temp)
  }
  
  
  rm(placebo_full_dr1_temp, placebo_full_dr2_temp, placebo_full_pr1_temp, placebo_full_pr2_temp)
}

#bind each run together
placebo_full <- rbind(placebo_full_1, placebo_full_2, placebo_full_3, placebo_full_4, placebo_full_5, placebo_full_6)
#indicate that these are all actually white jurors
placebo_full$nonwhite <- 0
#remove temp dataframes
rm(placebo_full_1, placebo_full_2, placebo_full_3, placebo_full_4, placebo_full_5, placebo_full_6)

#save out the subset for each group and its median gap
doll_rd1_placebo_full <- placebo_full %>% filter(round==1 & type=="dollar")
doll_rd1_placebo_full_med <- median(doll_rd1_placebo_full$diff.diff)

doll_rd2_placebo_full <- placebo_full %>% filter(round==2 & type=="dollar")
doll_rd2_placebo_full_med <- median(doll_rd2_placebo_full$diff.diff)

rate_rd1_placebo_full <- placebo_full %>% filter(round==1 & type=="rating")
rate_rd1_placebo_full_med <- median(rate_rd1_placebo_full$diff.diff)

rate_rd2_placebo_full <- placebo_full %>% filter(round==2 & type=="rating")
rate_rd2_placebo_full_med <- median(rate_rd2_placebo_full$diff.diff)

#combine actual and placebo dissenter dataframes
dissent_influence <- rbind(nonwhite_rep_actual, placebo_full)

#simple t-test of gaps between real nonwhite and placebo dissenters:
t.test(dissent_influence$diff.diff~dissent_influence$nonwhite)

#simple regression model 
nonwhite_influence_nofixed <- lm(diff.diff~nonwhite, data=dissent_influence)
summary(nonwhite_influence_nofixed)

#regression model with robust SEs, clusters, scenario controls
nonwhite_influence_robust_simple <- lm_robust(diff.diff~nonwhite+scenario.x, data=dissent_influence, clusters=jurynum)
summary(nonwhite_influence_robust_simple)
summary(margins(nonwhite_influence_robust_simple))

#calculate marginal diff and plot
emmfit <- emmeans(nonwhite_influence_robust_simple, "nonwhite")
emmfit_dt <- summary(emmfit) 
emmfit_dt$nonwhite <- factor(emmfit_dt$nonwhite, labels=c("White Dissenters", "POC Dissenters"))
emm_plot <- ggplot(data=emmfit_dt, aes(y=emmean, x=nonwhite, ymin = lower.CL, ymax = upper.CL) ) +
  expand_limits(y=c(0,2)) + 
  geom_point(size=4, stat="identity") +
  geom_errorbar(aes(), width=.05) +
  ylab("Estimated Influence Gap") +
  xlab("") +
  theme_bw()+
  theme(axis.title = element_text(size = 15), axis.text.x = element_text(size = 15))
emm_plot


## add matching on scenario and distance from other jurors
  #calculate difference in opinion
dissent_influence$init_distance <- dissent_influence$mean.y - dissent_influence$mean.x
  #find matches 
mod1 <- matchit(formula = nonwhite~init_distance+scenario.x,
                data = dissent_influence,
                method = "cem",
                estimand = "ATT", 
                na.action="na.omit")
  #add matching weights to dataframe
dissent_influence$w <- mod1$weights 

#add matching weights in model of influence gaps
nonwhite_influence_robust_matched <- lm_robust(diff.diff~nonwhite+scenario.x, data=dissent_influence, clusters=jurynum, weights=w)
emmfit_b <- emmeans(nonwhite_influence_robust_matched, "nonwhite")
emmfit_b <- summary(emmfit_b) 
emmfit_b$nonwhite <- factor(emmfit_b$nonwhite, labels=c("White Dissenters", "POC Dissenters"))
emmfit <- bind_rows(emmfit_dt, emmfit_b, .id="id") %>%
  mutate("Model" = case_when(id==1~"Unadjusted",
                             id==2~"With Matching"))
# add matched differences to plot
emm_plot <- ggplot(data=emmfit, aes(y=emmean, x=(nonwhite=="POC Dissenters")+(.1*(Model=="With Matching")), ymin = lower.CL, ymax = upper.CL, color=Model) ) +
  expand_limits(y=c(0,2)) + 
  geom_point(size=4, stat="identity") +
  geom_errorbar(aes(), width=.05) +
  ylab("Estimated Influence Gap") +
  xlab("") +
  theme_bw()+
  theme(axis.title = element_text(size = 15), axis.text.x = element_text(size = 15)) + 
  scale_color_manual(values=c("black", "gray60")) + 
  scale_x_continuous(breaks=c(0,1), labels=c("White Dissenters", "POC Dissenters")) + 
  coord_cartesian(xlim=c(-.5, 1.5))
emm_plot
ggsave("Figures/emm_plot.pdf")

rm(doll_rd1_placebo_full, doll_rd2_placebo_full, dr1_placebo_full, dr1_placebo.nonwhite_full, dr1_placebo.white_full,
   dr2_placebo_full, dr2_placebo.nonwhite_full, dr2_placebo.white_full,
   pr1_placebo_full, pr1_placebo.nonwhite_full, pr1_placebo.white_full,
   pr2_placebo_full, pr2_placebo.nonwhite_full, pr2_placebo.white_full,
   rate_rd1_placebo_full, rate_rd2_placebo_full)

## Permutation tests (incl. appendix figure A10) ####

#simple permutation test:
actual <- diff(by(dissent_influence$diff.diff, dissent_influence$nonwhite, mean))
dist <- replicate(2000, diff(by(dissent_influence$diff.diff, sample(dissent_influence$nonwhite, length(dissent_influence$nonwhite), FALSE), mean)))
dist_plot <-as.data.frame(dist)
quantile(dist_plot$dist, probs = seq(0, 1, 0.01))

plot_permutation <- ggplot(data=dist_plot)+
  geom_density(aes(dist),color="gray",fill="gray",alpha=0.2)+
  labs(x="Distribution of Permutations",y="Density",
       title="All Rounds Combined, Lone Dissenters")+
  geom_vline(xintercept = actual, color="blue") +
  expand_limits(x=c(-1,1)) +
  theme_bw()
plot_permutation

ggsave(here("Figures/plot_permutation.pdf"))

rm(actual, dist, dist_plot, plot_permutation)

# another method for permutation test: coin package (reported)
independence_test(diff.diff~nonwhite, data=dissent_influence)


## 9. Figure 4: influence gap for lone dissenters, white and nonwhite #####
  #set up and run loop ####
i <- 1

#create frames to store results from each round/scale
placebo_bootstrap1 <- data.frame(matrix(ncol = 1, nrow =0))
count_vars <- c("diffdiff")
colnames(placebo_bootstrap1) <- count_vars

placebo_bootstrap1_round1 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap1_round1) <- count_vars
placebo_bootstrap1_round2 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap1_round2) <- count_vars

placebo_bootstrap1_dr1_1 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap1_round1) <- count_vars
placebo_bootstrap1_dr2_1 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap1_round2) <- count_vars

placebo_bootstrap1_pr1_1 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap1_round1) <- count_vars
placebo_bootstrap1_pr2_1 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap1_round2) <- count_vars


set.seed(9904)

# repeating 500 times, randomly assign a placebo nonwhite dissenter and calculate influence
for (i in 1:500) {
  
  # assign placebo dissenter for pre-deliberation rounds
  placebo_a_bootstrap1 <- within(pre_del.homog, {
    random_a_bootstrap1 <- block_ra(blocks = pre_del.homog$jurynum, prob_each = c(.83, .17))
  })
  
  # assign placebo for post-deliberation rounds
  placebo_b_bootstrap1 <- within(post_del.homog, {
    random_b_bootstrap1 <- block_ra(blocks = post_del.homog$jurynum, prob_each = c(.83, .17))
  })
  
  #Calculate differences between placebo/other jurors (as in single-run version above)
  #Dollar Scale; Dollar-first (rd 1):
  dr1_1_placebo_bootstrap <- placebo_b_bootstrap1 %>% 
    group_by(jurynum,random_b_bootstrap1,jury.eight,scenario) %>% 
    summarize(mean = mean(eight.doll,na.rm=T))
  dr1_1_placebo.white_bootstrap <- dr1_1_placebo_bootstrap %>% filter(random_b_bootstrap1==0)
  dr1_1_placebo.nonwhite_bootstrap <- dr1_1_placebo_bootstrap %>% filter(random_b_bootstrap1==1)
  dr1_1_placebo_bootstrap <- left_join(dr1_1_placebo.white_bootstrap,dr1_1_placebo.nonwhite_bootstrap,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  summary(dr1_1_placebo_bootstrap$diff) 
  dr1_1_placebo_bootstrap <- dr1_1_placebo_bootstrap %>% mutate(diff=abs(diff))
  dr1_1_placebo_bootstrap$diff[is.nan(dr1_1_placebo_bootstrap$diff)]<-NA
  dr1_1_placebo_bootstrap$mean.x[is.nan(dr1_1_placebo_bootstrap$mean.x)]<-NA
  dr1_1_placebo_bootstrap$mean.y[is.nan(dr1_1_placebo_bootstrap$mean.y)]<-NA
  
  #label dissenters
  dr1_1_placebo_bootstrap$num_disagree_bootstrap <- 0
  dr1_1_placebo_bootstrap$num_disagree_bootstrap[dr1_1_placebo_bootstrap$diff>=1.5] <- 1

  
  #Dollar Scale; Rating-first (rd 2):
  dr2_1_placebo_bootstrap <- placebo_a_bootstrap1 %>% group_by(jurynum,random_a_bootstrap1,jury.eight,scenario) %>% 
    summarize(mean = mean(eight.doll,na.rm=T))
  dr2_1_placebo.white_bootstrap <- dr2_1_placebo_bootstrap %>% filter(random_a_bootstrap1==0)
  dr2_1_placebo.nonwhite_bootstrap <- dr2_1_placebo_bootstrap %>% filter(random_a_bootstrap1==1)
  dr2_1_placebo_bootstrap <- left_join(dr2_1_placebo.white_bootstrap,dr2_1_placebo.nonwhite_bootstrap,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  summary(dr2_1_placebo_bootstrap$diff) 
  dr2_1_placebo_bootstrap <- dr2_1_placebo_bootstrap %>% mutate(diff=abs(diff))
  dr2_1_placebo_bootstrap$diff[is.nan(dr2_1_placebo_bootstrap$diff)]<-NA
  dr2_1_placebo_bootstrap$mean.x[is.nan(dr2_1_placebo_bootstrap$mean.x)]<-NA
  dr2_1_placebo_bootstrap$mean.y[is.nan(dr2_1_placebo_bootstrap$mean.y)]<-NA
  summary(dr2_1_placebo_bootstrap$diff) 
  
  dr2_1_placebo_bootstrap$num_disagree <- 0
  dr2_1_placebo_bootstrap$num_disagree[dr2_1_placebo_bootstrap$diff>=1.5] <- 1
  table(dr2_1_placebo_bootstrap$num_disagree)  
  
  #Punishment Scale; Rating-first (rd 1):
  pr1_1_placebo_bootstrap <- placebo_a_bootstrap1 %>% group_by(jurynum,random_a_bootstrap1,jryscale,scenario) %>% 
    summarize(mean = mean(iscale,na.rm=T))
  pr1_1_placebo.white_bootstrap <- pr1_1_placebo_bootstrap %>% filter(random_a_bootstrap1==0)
  pr1_1_placebo.nonwhite_bootstrap <- pr1_1_placebo_bootstrap %>% filter(random_a_bootstrap1==1)
  pr1_1_placebo_bootstrap <- left_join(pr1_1_placebo.white_bootstrap,pr1_1_placebo.nonwhite_bootstrap,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  summary(pr1_1_placebo_bootstrap$diff) 
  pr1_1_placebo_bootstrap <- pr1_1_placebo_bootstrap %>% mutate(diff=abs(diff))
  pr1_1_placebo_bootstrap$diff[is.nan(pr1_1_placebo_bootstrap$diff)]<-NA
  pr1_1_placebo_bootstrap$mean.x[is.nan(pr1_1_placebo_bootstrap$mean.x)]<-NA
  pr1_1_placebo_bootstrap$mean.y[is.nan(pr1_1_placebo_bootstrap$mean.y)]<-NA
  summary(pr1_1_placebo_bootstrap$diff) 
  
  pr1_1_placebo_bootstrap$num_disagree <- 0
  pr1_1_placebo_bootstrap$num_disagree[pr1_1_placebo_bootstrap$diff>=1.5] <- 1
  table(pr1_1_placebo_bootstrap$num_disagree)  
  
  #Punishment Scale; Dollar-first (rd 2):
  pr2_1_placebo_bootstrap <- placebo_b_bootstrap1 %>% group_by(jurynum,random_b_bootstrap1,jryscale,scenario) %>% 
    summarize(mean = mean(iscale,na.rm=T))
  pr2_1_placebo.white_bootstrap <- pr2_1_placebo_bootstrap %>% filter(random_b_bootstrap1==0)
  pr2_1_placebo.nonwhite_bootstrap <- pr2_1_placebo_bootstrap %>% filter(random_b_bootstrap1==1)
  pr2_1_placebo_bootstrap <- left_join(pr2_1_placebo.white_bootstrap,pr2_1_placebo.nonwhite_bootstrap,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  summary(pr2_1_placebo_bootstrap$diff) 
  pr2_1_placebo_bootstrap <- pr2_1_placebo_bootstrap %>% mutate(diff=abs(diff))
  pr2_1_placebo_bootstrap$diff[is.nan(pr2_1_placebo_bootstrap$diff)]<-NA
  pr2_1_placebo_bootstrap$mean.x[is.nan(pr2_1_placebo_bootstrap$mean.x)]<-NA
  pr2_1_placebo_bootstrap$mean.y[is.nan(pr2_1_placebo_bootstrap$mean.y)]<-NA
  summary(pr2_1_placebo_bootstrap$diff) 
  
  pr2_1_placebo_bootstrap$num_disagree <- 0
  pr2_1_placebo_bootstrap$num_disagree[pr2_1_placebo_bootstrap$diff>=1.5] <- 1
  table(pr2_1_placebo_bootstrap$num_disagree)  
  
  sum_nonwhite1_placebo_bootstrap <- rbind(dr1_1_placebo_bootstrap, dr2_1_placebo_bootstrap, pr1_1_placebo_bootstrap, pr2_1_placebo_bootstrap)
  table(sum_nonwhite1_placebo_bootstrap$num_disagree)  
  
  
  ##Enforcing >1.5 Cutoffs 
  
  ##Dollar Scale; Dollar-first (rd 1); 1.5:
  dr1_1.a_placebo_bootstrap <- dr1_1_placebo_bootstrap %>%
    filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jury.eight.x),
           diff.nonwhite = abs(mean.y-jury.eight.x),
           diff.diff = diff.nonwhite-diff.white)
  doll_rd1_diff_placebo_bootstrap <- median(dr1_1.a_placebo_bootstrap$diff.diff, na.rm = T)
  
  
  ##Dollar Scale; Punishment-first (rd 2); 1.5:
  dr2_1.a_placebo_bootstrap <- dr2_1_placebo_bootstrap %>% 
    filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jury.eight.x),
           diff.nonwhite = abs(mean.y-jury.eight.x),
           diff.diff = diff.nonwhite-diff.white)
  doll_rd2_diff_placebo_bootstrap <- median(dr2_1.a_placebo_bootstrap$diff.diff)
  
  
  ##Punishment Scale; Punishment-first (rd 1); 1.5:
  pr1_1.a_placebo_bootstrap <- pr1_1_placebo_bootstrap %>% 
    filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jryscale.x),
           diff.nonwhite = abs(mean.y-jryscale.x),
           diff.diff = diff.nonwhite-diff.white)
  rating_rd1_diff_placebo_bootstrap <- median(pr1_1.a_placebo_bootstrap$diff.diff)
  
  
  ##Punishment Scale; Dollar-first (rd 2); 1.5:
  pr2_1.a_placebo_bootstrap <- pr2_1_placebo_bootstrap %>% 
    filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jryscale.x),
           diff.nonwhite = abs(mean.y-jryscale.x),
           diff.diff = diff.nonwhite-diff.white)
  rating_rd2_diff_placebo_bootstrap <- median(pr2_1.a_placebo_bootstrap$diff.diff)
  
  # save datasets for each round and calculate median gap
  jury_diffdiff_placebo_bootstrap <- rbind(dr1_1.a_placebo_bootstrap, dr2_1.a_placebo_bootstrap, pr1_1.a_placebo_bootstrap, pr2_1.a_placebo_bootstrap)
  diffdiff_med_placebo_bootstrap <- median(jury_diffdiff_placebo_bootstrap$diff.diff)
  
  jury_diffdiff_placebo_bootstrap_round1 <- rbind(dr1_1.a_placebo_bootstrap, pr1_1.a_placebo_bootstrap)
  diffdiff_med_placebo_bootstrap_round1 <- median(jury_diffdiff_placebo_bootstrap_round1$diff.diff)
  
  jury_diffdiff_placebo_bootstrap_round2 <- rbind(dr2_1.a_placebo_bootstrap, pr2_1.a_placebo_bootstrap)
  diffdiff_med_placebo_bootstrap_round2 <- median(jury_diffdiff_placebo_bootstrap_round2$diff.diff)
  
  # save differences in frames for each group
  temp.placebo <- data.frame(matrix(ncol = 1, nrow =1))
  colnames(temp.placebo) <- c("diff.diff")
  temp.placebo$diff.diff <- diffdiff_med_placebo_bootstrap
  
  temp.placebo_dr1_1 <- data.frame(matrix(ncol = 1, nrow = 1))
  colnames(temp.placebo_dr1_1) <- c("diff.diff")
  temp.placebo_dr1_1$diff.diff <- doll_rd1_diff_placebo_bootstrap
  
  temp.placebo_dr2_1 <- data.frame(matrix(ncol = 1, nrow = 1))
  colnames(temp.placebo_dr2_1) <- c("diff.diff")
  temp.placebo_dr2_1$diff.diff <- doll_rd2_diff_placebo_bootstrap
  
  temp.placebo_pr1_1 <- data.frame(matrix(ncol = 1, nrow = 1))
  colnames(temp.placebo_pr1_1) <- c("diff.diff")
  temp.placebo_pr1_1$diff.diff <- rating_rd1_diff_placebo_bootstrap
  
  temp.placebo_pr2_1 <- data.frame(matrix(ncol = 1, nrow = 1))
  colnames(temp.placebo_pr2_1) <- c("diff.diff")
  temp.placebo_pr2_1$diff.diff <- rating_rd2_diff_placebo_bootstrap
  
  temp.placebo_round1 <- data.frame(matrix(ncol = 1, nrow =1))
  colnames(temp.placebo_round1) <- c("diff.diff")
  temp.placebo_round1$diff.diff <- diffdiff_med_placebo_bootstrap_round1
  
  temp.placebo_round2 <- data.frame(matrix(ncol = 1, nrow =1))
  colnames(temp.placebo_round2) <- c("diff.diff")
  temp.placebo_round2$diff.diff <- diffdiff_med_placebo_bootstrap_round2
  
  # group new results together with previous runs
  placebo_bootstrap1 <- rbind(placebo_bootstrap1, temp.placebo)
  placebo_bootstrap1_dr1_1 <- rbind(placebo_bootstrap1_dr1_1, temp.placebo_dr1_1)
  placebo_bootstrap1_dr2_1 <- rbind(placebo_bootstrap1_dr2_1, temp.placebo_dr2_1)
  placebo_bootstrap1_pr1_1 <- rbind(placebo_bootstrap1_pr1_1, temp.placebo_pr1_1)
  placebo_bootstrap1_pr2_1 <- rbind(placebo_bootstrap1_pr2_1, temp.placebo_pr2_1)
  
  placebo_bootstrap1_round1 <- rbind(placebo_bootstrap1_round1, temp.placebo_round1)
  placebo_bootstrap1_round2 <- rbind(placebo_bootstrap1_round2, temp.placebo_round2)
  
}

  #plots (1 dissenter) ####

#create frame for storing results
placebo.data <- data.frame(matrix(ncol = 2, nrow =0))
placebo.data_round1 <- data.frame(matrix(ncol = 2, nrow =0))
placebo.data_round2 <- data.frame(matrix(ncol = 2, nrow =0))

#assign common variable names
main_vars <- c("jurynum", "diffdiff")
colnames(placebo.data) <- main_vars
colnames(placebo.data_round1) <- main_vars
colnames(placebo.data_round2) <- main_vars

# summarize simulation results
temp_dr1 <- dr1_placebo.a %>%
  group_by(jurynum) %>%
  summarize(diffdiff = median(diff.diff,na.rm=T))

temp_dr2 <- dr2_placebo.a %>%
  group_by(jurynum) %>%
  summarize(diffdiff = median(diff.diff,na.rm=T))

temp_pr1 <- pr1_placebo.a %>%
  group_by(jurynum) %>%
  summarize(diffdiff = median(diff.diff,na.rm=T))

temp_pr2 <- pr2_placebo.a %>%
  group_by(jurynum) %>%
  summarize(diffdiff = median(diff.diff,na.rm=T))

#combine group data
placebo.data <- rbind(placebo.data, temp_dr1)
placebo.data <- rbind(placebo.data, temp_dr2)
placebo.data <- rbind(placebo.data, temp_pr1)
placebo.data <- rbind(placebo.data, temp_pr2)

#calculate overall median
placebo_diffdiff_med <- median(placebo.data$diffdiff)

#combine real (non-placebo) dissenter groups and calculate differences
jury_diffdiff <- rbind(dr1.a, dr2.a, pr1.a, pr2.a)
diffdiff_med <- median(jury_diffdiff$diff.diff)
diffdiff_mean <- mean(jury_diffdiff$diff.diff)

#store differences by round/scale
doll_rd1_diff <- median(dr1.a$diff.diff, na.rm = T)
doll_rd2_diff <- median(dr2.a$diff.diff)
rating_rd1_diff <- median(pr1.a$diff.diff)
rating_rd2_diff <- median(pr2.a$diff.diff)

#plot distribution of gaps
plot_diffdiff_placebo_bootstrap<-ggplot(data=placebo_bootstrap1)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="All Rounds Combined, Lone Dissenters")+
  geom_vline(xintercept = diffdiff_med, color="blue") +
  expand_limits(x=c(-0.5,2)) +
  theme_bw()
plot_diffdiff_placebo_bootstrap

#view percentile gap sizes
quantile(placebo_bootstrap1$diff.diff, probs = seq(0, 1, 0.01))
#calculate rank
placebo_bootstrap1$perc_rank <- rank(placebo_bootstrap1$diff.diff)/length(placebo_bootstrap1$diff.diff)


# create datasets and figures for by-round plots 
jury_diffdiff_round1 <- rbind(dr1.a, pr1.a)
jury_diffdiff_round2 <- rbind(dr2.a, pr2.a)

diffdiff_med_round1 <- median(jury_diffdiff_round1$diff.diff)
diffdiff_mean_round1 <- mean(jury_diffdiff_round1$diff.diff)
diffdiff_sd_round1 <- sd(jury_diffdiff_round1$diff.diff)
diffdiff_se_round1 <- 1.253*(diffdiff_sd_round1/sqrt(length(jury_diffdiff_round1$diff.diff)))

diffdiff_med_round2 <- median(jury_diffdiff_round2$diff.diff)
diffdiff_mean_round2 <- mean(jury_diffdiff_round2$diff.diff)
diffdiff_sd_round2 <- sd(jury_diffdiff_round2$diff.diff)
diffdiff_se_round2 <- 1.253*(diffdiff_sd_round2/sqrt(length(jury_diffdiff_round2$diff.diff)))


## REPEAT for each round

plot_diffdiff_placebo_bootstrap_round1<-ggplot(data=placebo_bootstrap1_round1)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Round 1, Lone Dissenters")+
  geom_vline(xintercept = diffdiff_med_round1, color="blue") +
  expand_limits(x=c(-0.5,2.5)) +
  theme_bw()
plot_diffdiff_placebo_bootstrap_round1

#calculate percentile
quantile(placebo_bootstrap1_round1$diff.diff, probs = seq(0, 1, 0.01)) 
##Actual median 1.9: 98th percentile


plot_diffdiff_placebo_bootstrap_round2<-ggplot(data=placebo_bootstrap1_round2)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Round 2, Lone Dissenters")+
  geom_vline(xintercept = diffdiff_med_round2, color="blue") +
  expand_limits(x=c(-0.5,2.5)) +
  theme_bw()
plot_diffdiff_placebo_bootstrap_round2

#calculate median
quantile(placebo_bootstrap1_round2$diff.diff, probs = seq(0, 1, 0.01)) 
##Actual median 1.8: 91st percentile


# REPEAT for each round/scale combination

plot_diffdiff_placebo_bootstrap_dr1_1<-ggplot(data=placebo_bootstrap1_dr1_1)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Dollars, Round 1, Lone Dissenters")+
  geom_vline(xintercept = doll_rd1_diff, color="blue") +
  expand_limits(x=c(-0.5,3)) +
  theme_bw()
plot_diffdiff_placebo_bootstrap_dr1_1

quantile(placebo_bootstrap1_dr1_1$diff.diff, probs = seq(0, 1, 0.01)) 
##Actual median: 2.1, 96th percentile


plot_diffdiff_placebo_bootstrap_dr2_1<-ggplot(data=placebo_bootstrap1_dr2_1)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Dollars, Round 2, Lone Dissenters")+
  geom_vline(xintercept = doll_rd2_diff, color="blue") +
  expand_limits(x=c(-0.5,3)) +
  theme_bw()
plot_diffdiff_placebo_bootstrap_dr2_1

quantile(placebo_bootstrap1_dr2_1$diff.diff, probs = seq(0, 1, 0.01)) 
#Actual median: 2.4, 99th percentile


plot_diffdiff_placebo_bootstrap_pr1_1<-ggplot(data=placebo_bootstrap1_pr1_1)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Ratings, Round 1, Lone Dissenters")+
  geom_vline(xintercept = rating_rd1_diff, color="blue") +
  expand_limits(x=c(-0.5,3)) +
  theme_bw()
plot_diffdiff_placebo_bootstrap_pr1_1

quantile(placebo_bootstrap1_pr1_1$diff.diff, probs = seq(0, 1, 0.01)) 
#Actual median: 61-97th percentile

plot_diffdiff_placebo_bootstrap_pr2_1<-ggplot(data=placebo_bootstrap1_pr2_1)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Ratings, Round 2, Lone Dissenters")+
  geom_vline(xintercept = rating_rd2_diff, color="blue") +
  expand_limits(x=c(-0.5,3)) +
  theme_bw()
plot_diffdiff_placebo_bootstrap_pr2_1

quantile(placebo_bootstrap1_pr2_1$diff.diff, probs = seq(0, 1, 0.01)) 
#Actual median, 1.6: 37th-53rd percentile



  #create real data (2 dissenters) ####

# Subsets of 2-nonwhite juries:
pre_del.hetr2 <- pre_del %>% filter(totnwhite.inc==2)
post_del.hetr2 <- post_del %>% filter(totnwhite.inc==2)

#Calculating gaps for each subset

#Dollar Scale; Dollar-first (rd 1):
dr1_2 <- post_del.hetr2 %>% group_by(jurynum,ind_nonwhite,jury.eight,scenario) %>% 
  summarize(mean = mean(eight.doll,na.rm=T))
dr1_2.white <- dr1_2 %>% filter(ind_nonwhite==0)
dr1_2.nonwhite <- dr1_2 %>% filter(ind_nonwhite==1)
dr1_2 <- left_join(dr1_2.white,dr1_2.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
summary(dr1_2$diff) 
dr1_2 <- dr1_2 %>% mutate(diff=abs(diff))
dr1_2$diff[is.nan(dr1_2$diff)]<-NA
dr1_2$mean.x[is.nan(dr1_2$mean.x)]<-NA
dr1_2$mean.y[is.nan(dr1_2$mean.y)]<-NA

dr1_2$num_disagree <- 0
dr1_2$num_disagree[dr1_2$diff>=1.5] <- 1

#Dollar Scale; Rating-first (rd 2):
dr2_2 <- pre_del.hetr2 %>% group_by(jurynum,ind_nonwhite,jury.eight,scenario) %>% 
  summarize(mean = mean(eight.doll,na.rm=T))
dr2_2.white <- dr2_2 %>% filter(ind_nonwhite==0)
dr2_2.nonwhite <- dr2_2 %>% filter(ind_nonwhite==1)
dr2_2 <- left_join(dr2_2.white,dr2_2.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
dr2_2 <- dr2_2 %>% mutate(diff=abs(diff))
dr2_2$diff[is.nan(dr2_2$diff)]<-NA
dr2_2$mean.x[is.nan(dr2_2$mean.x)]<-NA
dr2_2$mean.y[is.nan(dr2_2$mean.y)]<-NA

dr2_2$num_disagree <- 0
dr2_2$num_disagree[dr2_2$diff>=1.5] <- 1

#Rating Scale; Rating-first (rd 1):
pr1_2 <- pre_del.hetr2 %>% group_by(jurynum,ind_nonwhite,jryscale,scenario) %>% 
  summarize(mean = mean(iscale,na.rm=T))
pr1_2.white <- pr1_2 %>% filter(ind_nonwhite==0)
pr1_2.nonwhite <- pr1_2 %>% filter(ind_nonwhite==1)
pr1_2 <- left_join(pr1_2.white,pr1_2.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
pr1_2 <- pr1_2 %>% mutate(diff=abs(diff))
pr1_2$diff[is.nan(pr1_2$diff)]<-NA
pr1_2$mean.x[is.nan(pr1_2$mean.x)]<-NA
pr1_2$mean.y[is.nan(pr1_2$mean.y)]<-NA

pr1_2$num_disagree <- 0
pr1_2$num_disagree[pr1_2$diff>=1.5] <- 1

#Rating Scale; Rating-first (rd 2):
pr2_2 <- post_del.hetr2 %>% group_by(jurynum,ind_nonwhite,jryscale,scenario) %>% 
  summarize(mean = mean(iscale,na.rm=T))
pr2_2.white <- pr2_2 %>% filter(ind_nonwhite==0)
pr2_2.nonwhite <- pr2_2 %>% filter(ind_nonwhite==1)
pr2_2 <- left_join(pr2_2.white,pr2_2.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
summary(pr2_2$diff) 
pr2_2 <- pr2_2 %>% mutate(diff=abs(diff))
pr2_2$diff[is.nan(pr2_2$diff)]<-NA
pr2_2$mean.x[is.nan(pr2_2$mean.x)]<-NA
pr2_2$mean.y[is.nan(pr2_2$mean.y)]<-NA

pr2_2$num_disagree <- 0
pr2_2$num_disagree[pr2_2$diff>=1.5] <- 1
table(pr2_2$num_disagree)  

sum_nonwhite2 <- rbind(dr1_2, dr2_2, pr1_2, pr2_2)

##Enforcing Cutoffs for Dissent:

##Dollar Scale; Dollar-first (rd 1); 1.5:
dr1_2.a <- dr1_2 %>%
  filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jury.eight.x),
         diff.nonwhite = abs(mean.y-jury.eight.x),
         diff.diff = diff.nonwhite-diff.white)
doll_rd1_diff2 <- median(dr1_2.a$diff.diff, na.rm = T)
#t.test(dr1.a$diff.white,dr1.a$diff.nonwhite) #Diff=0.7, p=0.08

##Dollar Scale; Rating-first (rd 2); 1.5:
dr2_2.a <- dr2_2 %>% 
  filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jury.eight.x),
         diff.nonwhite = abs(mean.y-jury.eight.x),
         diff.diff = diff.nonwhite-diff.white)
doll_rd2_diff2 <- median(dr2_2.a$diff.diff)
#t.test(dr2.a$diff.white,dr2.a$diff.nonwhite) #Diff=1.7 p=0.00

##Rating Scale; Rating-first (rd 1); 1.5:
pr1_2.a <- pr1_2 %>% 
  filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jryscale.x),
         diff.nonwhite = abs(mean.y-jryscale.x),
         diff.diff = diff.nonwhite-diff.white)
rating_rd1_diff2 <- median(pr1_2.a$diff.diff)

##Rating Scale; Dollar-first (rd 2); 1.5:
pr2_2.a <- pr2_2 %>% 
  filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jryscale.x),
         diff.nonwhite = abs(mean.y-jryscale.x),
         diff.diff = diff.nonwhite-diff.white)
rating_rd2_diff2 <- median(pr2_2.a$diff.diff)


# Combine data and summarize differences 
jury2_diffdiff <- rbind(dr1_2.a, dr2_2.a, pr1_2.a, pr2_2.a)
diffdiff_med2 <- median(jury2_diffdiff$diff.diff)
summary(jury2_diffdiff$diff.diff)
t.test(jury2_diffdiff$diff.nonwhite, jury2_diffdiff$diff.white)

#store data and results 
jury2_diffdiff_round1 <- rbind(dr1_2.a, pr1_2.a)
diffdiff_med2_round1 <- median(jury2_diffdiff_round1$diff.diff)

jury2_diffdiff_round2 <- rbind(dr2_2.a, pr2_2.a)
diffdiff_med2_round2 <- median(jury2_diffdiff_round2$diff.diff)

#######################################################################.
  ##Placebo setup: 2 nonwhite jurors (plus appendix figure A9) ####

set.seed(23214)

# Randomly assign placebo dissenters 
placebo_a <- within(pre_del.homog, {
  random_a <- block_ra(blocks = pre_del.homog$jurynum, prob_each = c(.67, .33))
})

placebo_b <- within(post_del.homog, {
  random_b <- block_ra(blocks = post_del.homog$jurynum, prob_each = c(.67, .33))
})

# Calculating placebo differences for each round 
#Dollar Scale; Dollar-first (rd 1):
dr1_2_placebo <- placebo_b %>% group_by(jurynum,random_b,jury.eight,scenario) %>% 
  summarize(mean = mean(eight.doll,na.rm=T))
dr1_2_placebo.white <- dr1_2_placebo %>% filter(random_b==0)
dr1_2_placebo.nonwhite <- dr1_2_placebo %>% filter(random_b==1)
dr1_2_placebo <- left_join(dr1_2_placebo.white,dr1_2_placebo.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
dr1_2_placebo <- dr1_2_placebo %>% mutate(diff=abs(diff))
dr1_2_placebo$diff[is.nan(dr1_2_placebo$diff)]<-NA
dr1_2_placebo$mean.x[is.nan(dr1_2_placebo$mean.x)]<-NA
dr1_2_placebo$mean.y[is.nan(dr1_2_placebo$mean.y)]<-NA

dr1_2_placebo$num_disagree <- 0
dr1_2_placebo$num_disagree[dr1_2_placebo$diff>=1.5] <- 1

#Dollar Scale; Rating-first (rd 2):
dr2_2_placebo <- placebo_a %>% group_by(jurynum,random_a,jury.eight,scenario) %>% 
  summarize(mean = mean(eight.doll,na.rm=T))
dr2_2_placebo.white <- dr2_2_placebo %>% filter(random_a==0)
dr2_2_placebo.nonwhite <- dr2_2_placebo %>% filter(random_a==1)
dr2_2_placebo <- left_join(dr2_2_placebo.white,dr2_2_placebo.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
dr2_2_placebo <- dr2_2_placebo %>% mutate(diff=abs(diff))
dr2_2_placebo$diff[is.nan(dr2_2_placebo$diff)]<-NA
dr2_2_placebo$mean.x[is.nan(dr2_2_placebo$mean.x)]<-NA
dr2_2_placebo$mean.y[is.nan(dr2_2_placebo$mean.y)]<-NA

dr2_2_placebo$num_disagree <- 0
dr2_2_placebo$num_disagree[dr2_2_placebo$diff>=1.5] <- 1

#Rating Scale; Rating-first (rd 1):
pr1_2_placebo <- placebo_a %>% group_by(jurynum,random_a,jryscale,scenario) %>% 
  summarize(mean = mean(iscale,na.rm=T))
pr1_2_placebo.white <- pr1_2_placebo %>% filter(random_a==0)
pr1_2_placebo.nonwhite <- pr1_2_placebo %>% filter(random_a==1)
pr1_2_placebo <- left_join(pr1_2_placebo.white,pr1_2_placebo.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
pr1_2_placebo <- pr1_2_placebo %>% mutate(diff=abs(diff))
pr1_2_placebo$diff[is.nan(pr1_2_placebo$diff)]<-NA
pr1_2_placebo$mean.x[is.nan(pr1_2_placebo$mean.x)]<-NA
pr1_2_placebo$mean.y[is.nan(pr1_2_placebo$mean.y)]<-NA

pr1_2_placebo$num_disagree <- 0
pr1_2_placebo$num_disagree[pr1_2_placebo$diff>=1.5] <- 1

#Rating Scale; Rating-first (rd 2):
pr2_2_placebo <- placebo_b %>% group_by(jurynum,random_b,jryscale,scenario) %>% 
  summarize(mean = mean(iscale,na.rm=T))
pr2_2_placebo.white <- pr2_2_placebo %>% filter(random_b==0)
pr2_2_placebo.nonwhite <- pr2_2_placebo %>% filter(random_b==1)
pr2_2_placebo <- left_join(pr2_2_placebo.white,pr2_2_placebo.nonwhite,by="jurynum") %>% mutate(diff = mean.x-mean.y)
summary(pr2_2_placebo$diff) 
pr2_2_placebo <- pr2_2_placebo %>% mutate(diff=abs(diff))
pr2_2_placebo$diff[is.nan(pr2_2_placebo$diff)]<-NA
pr2_2_placebo$mean.x[is.nan(pr2_2_placebo$mean.x)]<-NA
pr2_2_placebo$mean.y[is.nan(pr2_2_placebo$mean.y)]<-NA

pr2_2_placebo$num_disagree <- 0
pr2_2_placebo$num_disagree[pr2_2_placebo$diff>=1.5] <- 1

sum_nonwhite2_placebo <- rbind(dr1_2_placebo, dr2_2_placebo, pr1_2_placebo, pr2_2_placebo)


##Enforcing Cutoffs for Dissenter status
## (and plot differences)


##Dollar Scale; Dollar-first (rd 1); 1.5:
dr1_2.a_placebo <- dr1_2_placebo %>%
  filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jury.eight.x),
         diff.nonwhite = abs(mean.y-jury.eight.x),
         diff.diff = diff.nonwhite-diff.white)
doll_rd1_diff2_placebo <- median(dr1_2.a_placebo$diff.diff, na.rm = T)


##Dollar Scale; Punishment-first (rd 2); 1.5:
dr2_2.a_placebo <- dr2_2_placebo %>% 
  filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jury.eight.x),
         diff.nonwhite = abs(mean.y-jury.eight.x),
         diff.diff = diff.nonwhite-diff.white)
doll_rd2_diff2_placebo <- median(dr2_2.a_placebo$diff.diff)


##Punishment Scale; Punishment-first (rd 1); 1.5:
pr1_2.a_placebo <- pr1_2_placebo %>% 
  filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jryscale.x),
         diff.nonwhite = abs(mean.y-jryscale.x),
         diff.diff = diff.nonwhite-diff.white)
rating_rd1_diff2_placebo <- median(pr1_2.a_placebo$diff.diff)


plotpr1_2_placebo<-ggplot(data=pr1_2.a_placebo)+
  geom_density(aes(diff.white),color="red",fill="red",alpha=0.2)+
  geom_density(aes(diff.nonwhite),color="blue",fill="blue",alpha=0.2)+
  labs(x="Difference Between Group Mean and Jury Decision",y="Density",
       title="Rating; 1.5 Cutoff; Round 1 (Placebo)")+
  theme_bw()
plotpr1_2_placebo

##Punishment Scale; Dollar-first (rd 2); 1.5:
pr2_2.a_placebo <- pr2_2_placebo %>% 
  filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
  mutate(diff.white = abs(mean.x-jryscale.x),
         diff.nonwhite = abs(mean.y-jryscale.x),
         diff.diff = diff.nonwhite-diff.white)
rating_rd2_diff2_placebo <- median(pr2_2.a_placebo$diff.diff)

#store data and differences
jury2_diffdiff_placebo <- rbind(dr1_2.a_placebo, dr2_2.a_placebo, pr1_2.a_placebo, pr2_2.a_placebo)
diffdiff_med2_placebo <- median(jury2_diffdiff_placebo$diff.diff)

jury2_diffdiff_placebo_round1 <- rbind(dr1_2.a_placebo, pr1_2.a_placebo)
diffdiff_med2_placebo_round1 <- median(jury2_diffdiff_placebo_round1$diff.diff)

jury2_diffdiff_placebo_round2 <- rbind(dr2_2.a_placebo, pr2_2.a_placebo)
diffdiff_med2_placebo_round2 <- median(jury2_diffdiff_placebo_round2$diff.diff)



##Placebo test: 2 nonwhite jurors ####
i <- 1

#create frames to store results
placebo_bootstrap2 <- data.frame(matrix(ncol = 1, nrow =0))
count_vars <- c("diffdiff")
colnames(placebo_bootstrap2) <- count_vars

placebo_bootstrap2_round1 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap2_round1) <- count_vars

placebo_bootstrap2_round2 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap2_round2) <- count_vars

placebo_bootstrap2_dr1_2 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap2_round1) <- count_vars
placebo_bootstrap2_dr2_2 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap2_round2) <- count_vars

placebo_bootstrap2_pr1_2 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap2_round1) <- count_vars
placebo_bootstrap2_pr2_2 <- data.frame(matrix(ncol = 1, nrow =0))
colnames(placebo_bootstrap2_round2) <- count_vars

#run bootstrap (2 dissenters) 

set.seed(23214)

for (i in 1:500) {
  
  #randomly assign dissenters 
  placebo_a_bootstrap <- within(pre_del.homog, {
    random_a_bootstrap <- block_ra(blocks = pre_del.homog$jurynum, prob_each = c(.67, .33))
  })
  
  placebo_b_bootstrap <- within(post_del.homog, {
    random_b_bootstrap <- block_ra(blocks = post_del.homog$jurynum, prob_each = c(.67, .33))
  })
  
  
  #Calculating differences b/w dissenters and other jurors
  #Dollar Scale; Dollar-first (rd 1):
  dr1_2_placebo_bootstrap <- placebo_b_bootstrap %>% group_by(jurynum,random_b_bootstrap,jury.eight,scenario) %>% 
    summarize(mean = mean(eight.doll,na.rm=T))
  dr1_2_placebo.white_bootstrap <- dr1_2_placebo_bootstrap %>% filter(random_b_bootstrap==0)
  dr1_2_placebo.nonwhite_bootstrap <- dr1_2_placebo_bootstrap %>% filter(random_b_bootstrap==1)
  dr1_2_placebo_bootstrap <- left_join(dr1_2_placebo.white_bootstrap,dr1_2_placebo.nonwhite_bootstrap,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  summary(dr1_2_placebo_bootstrap$diff) 
  dr1_2_placebo_bootstrap <- dr1_2_placebo_bootstrap %>% mutate(diff=abs(diff))
  dr1_2_placebo_bootstrap$diff[is.nan(dr1_2_placebo_bootstrap$diff)]<-NA
  dr1_2_placebo_bootstrap$mean.x[is.nan(dr1_2_placebo_bootstrap$mean.x)]<-NA
  dr1_2_placebo_bootstrap$mean.y[is.nan(dr1_2_placebo_bootstrap$mean.y)]<-NA
  summary(dr1_2_placebo_bootstrap$diff) 
  
  dr1_2_placebo_bootstrap$num_disagree_bootstrap <- 0
  dr1_2_placebo_bootstrap$num_disagree_bootstrap[dr1_2_placebo_bootstrap$diff>=1.5] <- 1
  table(dr1_2_placebo_bootstrap$num_disagree)  
  
  #Dollar Scale; Rating-first (rd 2):
  dr2_2_placebo_bootstrap <- placebo_a_bootstrap %>% group_by(jurynum,random_a_bootstrap,jury.eight,scenario) %>% 
    summarize(mean = mean(eight.doll,na.rm=T))
  dr2_2_placebo.white_bootstrap <- dr2_2_placebo_bootstrap %>% filter(random_a_bootstrap==0)
  dr2_2_placebo.nonwhite_bootstrap <- dr2_2_placebo_bootstrap %>% filter(random_a_bootstrap==1)
  dr2_2_placebo_bootstrap <- left_join(dr2_2_placebo.white_bootstrap,dr2_2_placebo.nonwhite_bootstrap,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  summary(dr2_2_placebo_bootstrap$diff) 
  dr2_2_placebo_bootstrap <- dr2_2_placebo_bootstrap %>% mutate(diff=abs(diff))
  dr2_2_placebo_bootstrap$diff[is.nan(dr2_2_placebo_bootstrap$diff)]<-NA
  dr2_2_placebo_bootstrap$mean.x[is.nan(dr2_2_placebo_bootstrap$mean.x)]<-NA
  dr2_2_placebo_bootstrap$mean.y[is.nan(dr2_2_placebo_bootstrap$mean.y)]<-NA
  summary(dr2_2_placebo_bootstrap$diff) 
  
  dr2_2_placebo_bootstrap$num_disagree <- 0
  dr2_2_placebo_bootstrap$num_disagree[dr2_2_placebo_bootstrap$diff>=1.5] <- 1
  table(dr2_2_placebo_bootstrap$num_disagree)  
  
  #Ratings Scale; Rating-first (rd 1):
  pr1_2_placebo_bootstrap <- placebo_a_bootstrap %>% group_by(jurynum,random_a_bootstrap,jryscale,scenario) %>% 
    summarize(mean = mean(iscale,na.rm=T))
  pr1_2_placebo.white_bootstrap <- pr1_2_placebo_bootstrap %>% filter(random_a_bootstrap==0)
  pr1_2_placebo.nonwhite_bootstrap <- pr1_2_placebo_bootstrap %>% filter(random_a_bootstrap==1)
  pr1_2_placebo_bootstrap <- left_join(pr1_2_placebo.white_bootstrap,pr1_2_placebo.nonwhite_bootstrap,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  summary(pr1_2_placebo_bootstrap$diff) 
  pr1_2_placebo_bootstrap <- pr1_2_placebo_bootstrap %>% mutate(diff=abs(diff))
  pr1_2_placebo_bootstrap$diff[is.nan(pr1_2_placebo_bootstrap$diff)]<-NA
  pr1_2_placebo_bootstrap$mean.x[is.nan(pr1_2_placebo_bootstrap$mean.x)]<-NA
  pr1_2_placebo_bootstrap$mean.y[is.nan(pr1_2_placebo_bootstrap$mean.y)]<-NA
  summary(pr1_2_placebo_bootstrap$diff) 
  
  pr1_2_placebo_bootstrap$num_disagree <- 0
  pr1_2_placebo_bootstrap$num_disagree[pr1_2_placebo_bootstrap$diff>=1.5] <- 1
  table(pr1_2_placebo_bootstrap$num_disagree)  
  
  #Ratings Scale; Dollar-first (rd 2):
  pr2_2_placebo_bootstrap <- placebo_b_bootstrap %>% group_by(jurynum,random_b_bootstrap,jryscale,scenario) %>% 
    summarize(mean = mean(iscale,na.rm=T))
  pr2_2_placebo.white_bootstrap <- pr2_2_placebo_bootstrap %>% filter(random_b_bootstrap==0)
  pr2_2_placebo.nonwhite_bootstrap <- pr2_2_placebo_bootstrap %>% filter(random_b_bootstrap==1)
  pr2_2_placebo_bootstrap <- left_join(pr2_2_placebo.white_bootstrap,pr2_2_placebo.nonwhite_bootstrap,by="jurynum") %>% mutate(diff = mean.x-mean.y)
  summary(pr2_2_placebo_bootstrap$diff) 
  pr2_2_placebo_bootstrap <- pr2_2_placebo_bootstrap %>% mutate(diff=abs(diff))
  pr2_2_placebo_bootstrap$diff[is.nan(pr2_2_placebo_bootstrap$diff)]<-NA
  pr2_2_placebo_bootstrap$mean.x[is.nan(pr2_2_placebo_bootstrap$mean.x)]<-NA
  pr2_2_placebo_bootstrap$mean.y[is.nan(pr2_2_placebo_bootstrap$mean.y)]<-NA
  summary(pr2_2_placebo_bootstrap$diff) 
  
  pr2_2_placebo_bootstrap$num_disagree <- 0
  pr2_2_placebo_bootstrap$num_disagree[pr2_2_placebo_bootstrap$diff>=1.5] <- 1
  table(pr2_2_placebo_bootstrap$num_disagree)  
  
  sum_nonwhite2_placebo_bootstrap <- rbind(dr1_2_placebo_bootstrap, dr2_2_placebo_bootstrap, pr1_2_placebo_bootstrap, pr2_2_placebo_bootstrap)
  table(sum_nonwhite2_placebo_bootstrap$num_disagree)  
  
  
  ##Enforcing Cutoffs for dissent
  
  ##Dollar Scale; Dollar-first (rd 1); 1.5:
  dr1_2.a_placebo_bootstrap <- dr1_2_placebo_bootstrap %>%
    filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jury.eight.x),
           diff.nonwhite = abs(mean.y-jury.eight.x),
           diff.diff = diff.nonwhite-diff.white)
  doll_rd1_diff2_placebo_bootstrap <- median(dr1_2.a_placebo_bootstrap$diff.diff, na.rm = T)
  
  
  ##Dollar Scale; Rating-first (rd 2); 1.5:
  dr2_2.a_placebo_bootstrap <- dr2_2_placebo_bootstrap %>% 
    filter(diff>=1.5 & is.na(jury.eight.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jury.eight.x),
           diff.nonwhite = abs(mean.y-jury.eight.x),
           diff.diff = diff.nonwhite-diff.white)
  doll_rd2_diff2_placebo_bootstrap <- median(dr2_2.a_placebo_bootstrap$diff.diff)
  
  
  ##Rating Scale; Rating-first (rd 1); 1.5:
  pr1_2.a_placebo_bootstrap <- pr1_2_placebo_bootstrap %>% 
    filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jryscale.x),
           diff.nonwhite = abs(mean.y-jryscale.x),
           diff.diff = diff.nonwhite-diff.white)
  rating_rd1_diff2_placebo_bootstrap <- median(pr1_2.a_placebo_bootstrap$diff.diff)
  
  
  ##Rating Scale; Dollar-first (rd 2); 1.5:
  pr2_2.a_placebo_bootstrap <- pr2_2_placebo_bootstrap %>% 
    filter(diff>=1.5 & is.na(jryscale.x)==F) %>%   
    mutate(diff.white = abs(mean.x-jryscale.x),
           diff.nonwhite = abs(mean.y-jryscale.x),
           diff.diff = diff.nonwhite-diff.white)
  rating_rd2_diff2_placebo_bootstrap <- median(pr2_2.a_placebo_bootstrap$diff.diff)
  
  # combine rounds 
  jury2_diffdiff_placebo_bootstrap <- rbind(dr1_2.a_placebo_bootstrap, dr2_2.a_placebo_bootstrap, pr1_2.a_placebo_bootstrap, pr2_2.a_placebo_bootstrap)
  diffdiff_med2_placebo_bootstrap <- median(jury2_diffdiff_placebo_bootstrap$diff.diff)
  
  jury2_diffdiff_placebo_bootstrap_round1 <- rbind(dr1_2.a_placebo_bootstrap, pr1_2.a_placebo_bootstrap)
  diffdiff_med2_placebo_bootstrap_round1 <- median(jury2_diffdiff_placebo_bootstrap_round1$diff.diff)
  
  jury2_diffdiff_placebo_bootstrap_round2 <- rbind(dr2_2.a_placebo_bootstrap, pr2_2.a_placebo_bootstrap)
  diffdiff_med2_placebo_bootstrap_round2 <- median(jury2_diffdiff_placebo_bootstrap_round2$diff.diff)
  
  #store results of this round and bind them with the results of the others
  temp.placebo <- data.frame(matrix(ncol = 1, nrow =1))
  colnames(temp.placebo) <- c("diff.diff")
  temp.placebo$diff.diff <- diffdiff_med2_placebo_bootstrap
  
  placebo_bootstrap2 <- rbind(placebo_bootstrap2, temp.placebo)
  
  temp.placebo_dr1_2 <- data.frame(matrix(ncol = 1, nrow = 1))
  colnames(temp.placebo_dr1_2) <- c("diff.diff")
  temp.placebo_dr1_2$diff.diff <- doll_rd1_diff2_placebo_bootstrap
  
  temp.placebo_dr2_2 <- data.frame(matrix(ncol = 1, nrow = 1))
  colnames(temp.placebo_dr2_2) <- c("diff.diff")
  temp.placebo_dr2_2$diff.diff <- doll_rd2_diff2_placebo_bootstrap
  
  temp.placebo_pr1_2 <- data.frame(matrix(ncol = 1, nrow = 1))
  colnames(temp.placebo_pr1_2) <- c("diff.diff")
  temp.placebo_pr1_2$diff.diff <- rating_rd1_diff2_placebo_bootstrap
  
  temp.placebo_pr2_2 <- data.frame(matrix(ncol = 1, nrow = 1))
  colnames(temp.placebo_pr2_2) <- c("diff.diff")
  temp.placebo_pr2_2$diff.diff <- rating_rd2_diff2_placebo_bootstrap
  
  placebo_bootstrap2_dr1_2 <- rbind(placebo_bootstrap2_dr1_2, temp.placebo_dr1_2)
  placebo_bootstrap2_dr2_2 <- rbind(placebo_bootstrap2_dr2_2, temp.placebo_dr2_2)
  placebo_bootstrap2_pr1_2 <- rbind(placebo_bootstrap2_pr1_2, temp.placebo_pr1_2)
  placebo_bootstrap2_pr2_2 <- rbind(placebo_bootstrap2_pr2_2, temp.placebo_pr2_2)
  
  
  temp.placebo_round1 <- data.frame(matrix(ncol = 1, nrow =1))
  colnames(temp.placebo_round1) <- c("diff.diff")
  temp.placebo_round1$diff.diff <- diffdiff_med2_placebo_bootstrap_round1
  
  placebo_bootstrap2_round1 <- rbind(placebo_bootstrap2_round1, temp.placebo_round1)
  
  temp.placebo_round2 <- data.frame(matrix(ncol = 1, nrow =1))
  colnames(temp.placebo_round2) <- c("diff.diff")
  temp.placebo_round2$diff.diff <- diffdiff_med2_placebo_bootstrap_round2
  
  placebo_bootstrap2_round2 <- rbind(placebo_bootstrap2_round2, temp.placebo_round2)
  
}

# grab bootstrap values and create figs ####
placebo_bootstrap2_med <- median(placebo_bootstrap2$diff.diff)
placebo_bootstrap2_sd <- sd(placebo_bootstrap2$diff.diff)
placebo_bootstrap2_se <- 1.253*(sd(placebo_bootstrap2$diff.diff)/sqrt(length(placebo_bootstrap2$diff.diff)))

placebo_bootstrap2_med_round1 <- median(placebo_bootstrap2_round1$diff.diff)
placebo_bootstrap2_sd_round1 <- sd(placebo_bootstrap2_round1$diff.diff)
placebo_bootstrap2_se_round1 <- 1.253*(sd(placebo_bootstrap2_round1$diff.diff)/sqrt(length(placebo_bootstrap2_round1$diff.diff)))

placebo_bootstrap2_med_round2 <- median(placebo_bootstrap2_round2$diff.diff)
placebo_bootstrap2_sd_round2 <- sd(placebo_bootstrap2_round2$diff.diff)
placebo_bootstrap2_se_round2 <- 1.253*(sd(placebo_bootstrap2_round2$diff.diff)/sqrt(length(placebo_bootstrap2_round2$diff.diff)))

# OVERALL difference plot:
plot_diffdiff2_placebo_bootstrap<-ggplot(data=placebo_bootstrap2)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="All Rounds Combined, Two Dissenters")+
  geom_vline(xintercept = diffdiff_med2, color="blue") +
  expand_limits(x=c(-0.5,2)) +
  theme_bw()
plot_diffdiff2_placebo_bootstrap

#calculate gap in percentiles
quantile(placebo_bootstrap2$diff.diff, probs = seq(0, 1, 0.01))
#Actual median: 1.5, 99th percentile

##Results for each round, combining DVs
plot_diffdiff2_placebo_bootstrap_round1<-ggplot(data=placebo_bootstrap2_round1)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Round 1, Two Dissenters")+
  geom_vline(xintercept = diffdiff_med2_round1, color="blue") +
  expand_limits(x=c(-0.5,2)) +
  theme_bw()
plot_diffdiff2_placebo_bootstrap_round1

quantile(placebo_bootstrap2_round1$diff.diff, probs = seq(0, 1, 0.01))
#Actual median: 1.5, 98th percentile


plot_diffdiff2_placebo_bootstrap_round2<-ggplot(data=placebo_bootstrap2_round2)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Round 2, Two Dissenters")+
  geom_vline(xintercept = diffdiff_med2_round2, color="blue") +
  expand_limits(x=c(-0.5,2)) +
  theme_bw()
plot_diffdiff2_placebo_bootstrap_round2

quantile(placebo_bootstrap2_round2$diff.diff, probs = seq(0, 1, 0.01))
#Actual median: 1.75, 99th percentile



##Now for each DV and round separately
plot_diffdiff2_placebo_bootstrap_dr1_2<-ggplot(data=placebo_bootstrap2_dr1_2)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Dollars, Round 1, Two Dissenters")+
  geom_vline(xintercept = doll_rd1_diff2, color="blue") +
  expand_limits(x=c(-0.5,3)) +
  theme_bw()
plot_diffdiff2_placebo_bootstrap_dr1_2

quantile(placebo_bootstrap2_dr1_2$diff.diff, probs = seq(0, 1, 0.01))
#Actual Nonwhite median: 1.5 (55th-79th percentile)


plot_diffdiff2_placebo_bootstrap_dr2_2<-ggplot(data=placebo_bootstrap2_dr2_2)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Dollars, Round 2, Two Dissenters")+
  geom_vline(xintercept = doll_rd2_diff2, color="blue") +
  expand_limits(x=c(-0.5,3)) +
  theme_bw()
plot_diffdiff2_placebo_bootstrap_dr2_2

quantile(placebo_bootstrap2_dr2_2$diff.diff, probs = seq(0, 1, 0.01))
#Actual Nonwhite median: 2 (99th+ percentile)


plot_diffdiff2_placebo_bootstrap_pr1_2<-ggplot(data=placebo_bootstrap2_pr1_2)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Ratings, Round 1, Two Dissenters")+
  geom_vline(xintercept = rating_rd1_diff2, color="blue") +
  expand_limits(x=c(-0.5,3)) +
  theme_bw()
plot_diffdiff2_placebo_bootstrap_pr1_2

quantile(placebo_bootstrap2_pr1_2$diff.diff, probs = seq(0, 1, 0.01))
#Actual Nonwhite median: 1 (62-71st percentile)



plot_diffdiff2_placebo_bootstrap_pr2_2<-ggplot(data=placebo_bootstrap2_pr2_2)+
  geom_density(aes(diff.diff),color="gray",fill="gray",alpha=0.2)+
  labs(x="Median Influence Gap",y="Density",
       title="Ratings, Round 2, Two Dissenters")+
  geom_vline(xintercept = rating_rd2_diff2, color="blue") +
  expand_limits(x=c(-0.5,3)) +
  theme_bw()
plot_diffdiff2_placebo_bootstrap_pr2_2

quantile(placebo_bootstrap2_pr2_2$diff.diff, probs = seq(0, 1, 0.01))
#Actual Nonwhite median: 1.33 (76.5th percentile)



##Create figures for paper ##

##All Rounds Aggregated
ggsave("Figures/placebo_figs_allrounds.pdf",arrangeGrob(plot_diffdiff_placebo_bootstrap, plot_diffdiff2_placebo_bootstrap, ncol=1,
                                                                       top = textGrob("",gp=gpar(fontsize=20,font=3)),
                                                                       bottom = textGrob("",gp=gpar(fontsize=8,font=3))),width=14,height=10)


##Round by Round (Appendix)
##Rounds separately
ggsave("Appendix Figures/placebo_figs_byround.pdf",arrangeGrob(plot_diffdiff_placebo_bootstrap_round1, plot_diffdiff2_placebo_bootstrap_round1, plot_diffdiff_placebo_bootstrap_round2, plot_diffdiff2_placebo_bootstrap_round2,
                                                                     ncol=2,
                                                                     top = textGrob("",gp=gpar(fontsize=20,font=3)),
                                                                     bottom = textGrob("",gp=gpar(fontsize=8,font=3))),width=14,height=10)


##Each DV separately
ggsave("Figures/placebo_figs_byround_bydv.pdf",arrangeGrob(plot_diffdiff_placebo_bootstrap_pr1_1, plot_diffdiff_placebo_bootstrap_pr2_1, plot_diffdiff_placebo_bootstrap_dr1_1, plot_diffdiff_placebo_bootstrap_dr2_1, plot_diffdiff2_placebo_bootstrap_pr1_2, plot_diffdiff2_placebo_bootstrap_pr2_2, plot_diffdiff2_placebo_bootstrap_dr1_2, plot_diffdiff2_placebo_bootstrap_dr2_2, ncol=2,
                                                                          top = textGrob("",gp=gpar(fontsize=20,font=3)),
                                                                          bottom = textGrob("",gp=gpar(fontsize=8,font=3))),width=14,height=10)




